Morphological variation of amazon and sailfin mollies

In this study, I am looking at various populations of amazon and sailfin mollies across their native/introduced range to assess morphological variation both within and among the species. I am using the Pickle fish collections for my samples.

All data

I will use all data first to see the trends, then filter for confirmed adult sizes (see “Filtered”).

Initial steps

Data collection

## Using libcurl 8.3.0 with Schannel

Checking normality

Checking normality

I will use Shapiro-wilke, histograms, and QQ plots to determine what traits are normal. These will only be performed on continuous variables, as discrete variables are not normal by nature.

Conclusions: literally all of them are NOT normal… will log transform them and run parametric tests (t-test, F-test, ANOVA). Not sure what to do with the discrete variables…[Logan said to use non-parametric OR a glmer with poisson distribution]

SW Test
shapiro.test(raw1$SL)
## 
##  Shapiro-Wilk normality test
## 
## data:  raw1$SL
## W = 0.9763, p-value = 2.763e-05
shapiro.test(raw1$BD)
## 
##  Shapiro-Wilk normality test
## 
## data:  raw1$BD
## W = 0.96287, p-value = 1.787e-07
shapiro.test(raw1$CPD)
## 
##  Shapiro-Wilk normality test
## 
## data:  raw1$CPD
## W = 0.96431, p-value = 2.908e-07
shapiro.test(raw1$CPL)
## 
##  Shapiro-Wilk normality test
## 
## data:  raw1$CPL
## W = 0.97288, p-value = 6.831e-06
shapiro.test(raw1$PreDL)
## 
##  Shapiro-Wilk normality test
## 
## data:  raw1$PreDL
## W = 0.97924, p-value = 9.997e-05
shapiro.test(raw1$DbL)
## 
##  Shapiro-Wilk normality test
## 
## data:  raw1$DbL
## W = 0.97697, p-value = 3.68e-05
shapiro.test(raw1$HL)
## 
##  Shapiro-Wilk normality test
## 
## data:  raw1$HL
## W = 0.94955, p-value = 3.057e-09
shapiro.test(raw1$HD)
## 
##  Shapiro-Wilk normality test
## 
## data:  raw1$HD
## W = 0.96761, p-value = 9.28e-07
shapiro.test(raw1$HW)
## 
##  Shapiro-Wilk normality test
## 
## data:  raw1$HW
## W = 0.97201, p-value = 4.833e-06
shapiro.test(raw1$SnL)
## 
##  Shapiro-Wilk normality test
## 
## data:  raw1$SnL
## W = 0.65042, p-value < 2.2e-16
shapiro.test(raw1$OL)
## 
##  Shapiro-Wilk normality test
## 
## data:  raw1$OL
## W = 0.98631, p-value = 0.003102
Histograms
hist(raw1$SL)

hist(raw1$BD)

hist(raw1$CPD)

hist(raw1$CPL)

hist(raw1$PreDL)

hist(raw1$DbL)

hist(raw1$HL)

hist(raw1$HD)

hist(raw1$HW)

hist(raw1$SnL)

hist(raw1$OL)

QQ plots
qqnorm(raw1$SL)
qqline(raw1$SL)

qqnorm(raw1$BD)
qqline(raw1$BD)

qqnorm(raw1$CPD)
qqline(raw1$CPD)

qqnorm(raw1$CPL)
qqline(raw1$CPL)

qqnorm(raw1$PreDL)
qqline(raw1$PreDL)

qqnorm(raw1$DbL)
qqline(raw1$DbL)

qqnorm(raw1$HL)
qqline(raw1$HL)

qqnorm(raw1$HD)
qqline(raw1$HD)

qqnorm(raw1$HW)
qqline(raw1$HW)

qqnorm(raw1$SnL)
qqline(raw1$SnL)

qqnorm(raw1$OL)
qqline(raw1$OL)

Log transformations

Since all of the continuous characters were not normal, I will log transform them, and proceed with the transformed data in all analyses.

raw1[paste0(names(raw1), '_log')] <- log(raw1[, 26:35], 10)
#for some reason did the log of the whole data frame instead of those specific columns, so I will just remove the unnecessary columns instead of spending time on what went wrong.

raw1 <- raw1[-c(37:60)]
raw1 <- raw1[-c(48)]

#double checking normality with a SW test and QQ plot

shapiro.test(raw1$SL_log)
## 
##  Shapiro-Wilk normality test
## 
## data:  raw1$SL_log
## W = 0.99271, p-value = 0.1056
shapiro.test(raw1$BD_log)
## 
##  Shapiro-Wilk normality test
## 
## data:  raw1$BD_log
## W = 0.82598, p-value < 2.2e-16
shapiro.test(raw1$CPD_log)
## 
##  Shapiro-Wilk normality test
## 
## data:  raw1$CPD_log
## W = 0.99285, p-value = 0.1144
shapiro.test(raw1$CPL_log)
## 
##  Shapiro-Wilk normality test
## 
## data:  raw1$CPL_log
## W = 0.99513, p-value = 0.3824
shapiro.test(raw1$PreDL_log)
## 
##  Shapiro-Wilk normality test
## 
## data:  raw1$PreDL_log
## W = 0.93554, p-value = 8.204e-11
shapiro.test(raw1$DbL_log)
## 
##  Shapiro-Wilk normality test
## 
## data:  raw1$DbL_log
## W = 0.98042, p-value = 0.0001711
shapiro.test(raw1$HL_log)
## 
##  Shapiro-Wilk normality test
## 
## data:  raw1$HL_log
## W = 0.99386, p-value = 0.1984
shapiro.test(raw1$HD_log)
## 
##  Shapiro-Wilk normality test
## 
## data:  raw1$HD_log
## W = 0.99661, p-value = 0.7102
shapiro.test(raw1$HW_log)
## 
##  Shapiro-Wilk normality test
## 
## data:  raw1$HW_log
## W = 0.99491, p-value = 0.3435
shapiro.test(raw1$SnL_log)
## 
##  Shapiro-Wilk normality test
## 
## data:  raw1$SnL_log
## W = 0.99366, p-value = 0.1788
shapiro.test(raw1$OL_log)
## 
##  Shapiro-Wilk normality test
## 
## data:  raw1$OL_log
## W = 0.99271, p-value = 0.1056
qqnorm(raw1$SL_log)
qqline(raw1$SL_log)

qqnorm(raw1$BD_log)
qqline(raw1$BD_log)

qqnorm(raw1$CPD_log)
qqline(raw1$CPD_log)

qqnorm(raw1$CPL_log)
qqline(raw1$CPL_log)

qqnorm(raw1$PreDL_log)
qqline(raw1$PreDL_log)

qqnorm(raw1$DbL_log)
qqline(raw1$DbL_log)

qqnorm(raw1$HL_log)
qqline(raw1$HL_log)

qqnorm(raw1$HD_log)
qqline(raw1$HD_log)

qqnorm(raw1$HW_log)
qqline(raw1$HW_log)

qqnorm(raw1$SnL_log)
qqline(raw1$SnL_log)

qqnorm(raw1$OL_log)
qqline(raw1$OL_log)

Correcting for body size (residuals)

SL correction

Since amazons are in general bigger than sailfin, we don’t want any results to be due to this difference in body size bias. Therefore, we will see what traits are influenced by body size (regressions) and correct for body size when necessary (absolute value of residuals). We can then use the residuals when comparing between species for traits that are influenced by body size, and raw/transformed data for traits that are not influenced by body size. I will also calculate standardized residuals to compare residuals across traits in later analyses.

Quick results summary: traits not influenced by body size are left & right pelvic, anal, scales above and below lateral line (except in mexicana), scales before dorsal fin and fluctuating asymmetry; all other traits influenced by body size.

Regressions

  • STEP ONE: plot trait vs body size
library(ggplot2)
library(ggpubr)



lat <- raw1[raw1$SPP == "p.latipinna",]


form <- raw1[raw1$SPP == "p.formosa",]


mex <- raw1[raw1$SPP == "p.mexicana",]



##### LAT #####
reg.lat.D <- lm(lat$D ~ lat$SL)
sd.lat.D <- rstandard(reg.lat.D)
reg.lat.D.plot <- ggplot(lat, aes(x = SL, y = D)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
reg.lat.D.plot
## `geom_smooth()` using formula = 'y ~ x'

reg.lat.P1 <- lm(lat$P1 ~ lat$SL)
sd.lat.P1 <- rstandard(reg.lat.P1)
reg.lat.P1.plot <- ggplot(lat, aes(x = SL, y = P1)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
reg.lat.P1.plot
## `geom_smooth()` using formula = 'y ~ x'

reg.lat.P2.L <- lm(lat$P2.L ~ lat$SL)
sd.lat.P2.L <- rstandard(reg.lat.P2.L)
reg.lat.P2.L.plot <- ggplot(lat, aes(x = SL, y = P2.L)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
reg.lat.P2.L.plot
## `geom_smooth()` using formula = 'y ~ x'

reg.lat.P2.R <- lm(lat$P2.R ~ lat$SL)
sd.lat.P2.R <- rstandard(reg.lat.P2.R)
reg.lat.P2.R.plot <- ggplot(lat, aes(x = SL, y = P2.R)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
reg.lat.P2.R.plot
## `geom_smooth()` using formula = 'y ~ x'

reg.lat.A <- lm(lat$A ~ lat$SL)
sd.lat.A <- rstandard(reg.lat.A)
reg.lat.A.plot <- ggplot(lat, aes(x = SL, y = A)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
reg.lat.A.plot
## `geom_smooth()` using formula = 'y ~ x'

reg.lat.P1.R <- lm(lat$P1.R ~ lat$SL)
sd.lat.P1.R <- rstandard(reg.lat.P1.R)
reg.lat.P1.R.plot <- ggplot(lat, aes(x = SL, y = P1.R)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
reg.lat.P1.R.plot
## `geom_smooth()` using formula = 'y ~ x'

reg.lat.LLSC <- lm(lat$LLSC ~ lat$SL)
sd.lat.LLSC <- rstandard(reg.lat.LLSC)
reg.lat.LLSC.plot <- ggplot(lat, aes(x = SL, y = LLSC)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
reg.lat.LLSC.plot
## `geom_smooth()` using formula = 'y ~ x'

reg.lat.SALL <- lm(lat$SALL ~ lat$SL)
sd.lat.SALL <- rstandard(reg.lat.SALL)
reg.lat.SALL.plot <- ggplot(lat, aes(x = SL, y = SALL)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
reg.lat.SALL.plot
## `geom_smooth()` using formula = 'y ~ x'

reg.lat.SBLL <- lm(lat$SBLL ~ lat$SL)
sd.lat.SBLL <- rstandard(reg.lat.SBLL)
reg.lat.SBLL.plot <- ggplot(lat, aes(x = SL, y = SBLL)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
reg.lat.SBLL.plot
## `geom_smooth()` using formula = 'y ~ x'

reg.lat.SBDF <- lm(lat$SBDF ~ lat$SL)
sd.lat.SBDF <- rstandard(reg.lat.SBDF)
reg.lat.SBDF.plot <- ggplot(lat, aes(x = SL, y = SBDF)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
reg.lat.SBDF.plot
## `geom_smooth()` using formula = 'y ~ x'

reg.lat.BD <- lm(lat$BD_log ~ lat$SL_log)
sd.lat.BD <- rstandard(reg.lat.BD)
reg.lat.BD.plot <- ggplot(lat, aes(x = SL_log, y = BD_log)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
reg.lat.BD.plot
## `geom_smooth()` using formula = 'y ~ x'

reg.lat.CPD <- lm(lat$CPD_log ~ lat$SL_log)
sd.lat.CPD <- rstandard(reg.lat.CPD)
reg.lat.CPD.plot <- ggplot(lat, aes(x = SL_log, y = CPD_log)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
reg.lat.CPD.plot
## `geom_smooth()` using formula = 'y ~ x'

reg.lat.CPL <- lm(lat$CPL_log ~ lat$SL_log)
sd.lat.CPL <- rstandard(reg.lat.CPL)
reg.lat.CPL.plot <- ggplot(lat, aes(x = SL_log, y = CPL_log)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
reg.lat.CPL.plot
## `geom_smooth()` using formula = 'y ~ x'

reg.lat.PreDL <- lm(lat$PreDL_log ~ lat$SL_log)
sd.lat.PreDL <- rstandard(reg.lat.PreDL)
reg.lat.PreDL.plot <- ggplot(lat, aes(x = SL_log, y = PreDL_log)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
reg.lat.PreDL.plot
## `geom_smooth()` using formula = 'y ~ x'

reg.lat.DbL <- lm(lat$DbL_log ~ lat$SL_log)
sd.lat.DbL <- rstandard(reg.lat.DbL)
reg.lat.DbL.plot <- ggplot(lat, aes(x = SL_log, y = DbL_log)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
reg.lat.DbL.plot
## `geom_smooth()` using formula = 'y ~ x'

reg.lat.HL <- lm(lat$HL_log ~ lat$SL_log)
sd.lat.HL <- rstandard(reg.lat.HL)
reg.lat.HL.plot <- ggplot(lat, aes(x = SL_log, y = HL_log)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
reg.lat.HL.plot
## `geom_smooth()` using formula = 'y ~ x'

reg.lat.HD <- lm(lat$HD_log ~ lat$SL_log)
sd.lat.HD <- rstandard(reg.lat.HD)
reg.lat.HD.plot <- ggplot(lat, aes(x = SL_log, y = HD_log)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
reg.lat.HD.plot
## `geom_smooth()` using formula = 'y ~ x'

reg.lat.HW <- lm(lat$HW_log ~ lat$SL_log)
sd.lat.HW <- rstandard(reg.lat.HW)
reg.lat.HW.plot <- ggplot(lat, aes(x = SL_log, y = HW_log)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
reg.lat.HW.plot
## `geom_smooth()` using formula = 'y ~ x'

reg.lat.SnL <- lm(lat$SnL_log ~ lat$SL_log)
sd.lat.SnL <- rstandard(reg.lat.SnL)
reg.lat.SnL.plot <- ggplot(lat, aes(x = SL_log, y = SnL_log)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
reg.lat.SnL.plot
## `geom_smooth()` using formula = 'y ~ x'

reg.lat.OL <- lm(lat$OL_log ~ lat$SL_log)
sd.lat.OL <- rstandard(reg.lat.OL)
reg.lat.OL.plot <- ggplot(lat, aes(x = SL_log, y = OL_log)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
reg.lat.OL.plot
## `geom_smooth()` using formula = 'y ~ x'

reg.lat.FLA <- lm(lat$FLA ~ lat$SL)
sd.lat.FLA <- rstandard(reg.lat.FLA)
reg.lat.FLA.plot <- ggplot(lat, aes(x = SL, y = FLA)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
reg.lat.FLA.plot
## `geom_smooth()` using formula = 'y ~ x'

##### FORM #####

reg.form.D <- lm(form$D ~ form$SL)
sd.form.D <- rstandard(reg.form.D)
reg.form.D.plot <- ggplot(form, aes(x =SL, y = D)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
reg.form.D.plot
## `geom_smooth()` using formula = 'y ~ x'

reg.form.P1 <- lm(form$P1 ~ form$SL)
sd.form.P1 <- rstandard(reg.form.P1)
reg.form.P1.plot <- ggplot(form, aes(x = SL, y = P1)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
reg.form.P1.plot
## `geom_smooth()` using formula = 'y ~ x'

reg.form.P2.L <- lm(form$P2.L ~ form$SL)
sd.form.P2.L <- rstandard(reg.form.P2.L)
reg.form.P2.L.plot <- ggplot(form, aes(x = SL, y = P2.L)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
reg.form.P2.L.plot
## `geom_smooth()` using formula = 'y ~ x'

reg.form.P2.R <- lm(form$P2.R ~ form$SL)
sd.form.P2.R <- rstandard(reg.form.P2.R)
reg.form.P2.R.plot <- ggplot(form, aes(x = SL, y = P2.R)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
reg.form.P2.R.plot
## `geom_smooth()` using formula = 'y ~ x'

reg.form.A <- lm(form$A ~ form$SL)
sd.form.A <- rstandard(reg.form.A)
reg.form.A.plot <- ggplot(form, aes(x = SL, y = A)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
reg.form.A.plot
## `geom_smooth()` using formula = 'y ~ x'

reg.form.P1.R <- lm(form$P1.R ~ form$SL)
sd.form.P1.R <- rstandard(reg.form.P1.R)
reg.form.P1.R.plot <- ggplot(form, aes(x = SL, y = P1.R)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
reg.form.P1.R.plot
## `geom_smooth()` using formula = 'y ~ x'

reg.form.LLSC <- lm(form$LLSC ~ form$SL)
sd.form.LLSC <- rstandard(reg.form.LLSC)
reg.form.LLSC.plot <- ggplot(form, aes(x = SL, y = LLSC)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
reg.form.LLSC.plot
## `geom_smooth()` using formula = 'y ~ x'

reg.form.SALL <- lm(form$SALL ~ form$SL)
sd.form.SALL <- rstandard(reg.form.SALL)
reg.form.SALL.plot <- ggplot(form, aes(x = SL, y = SALL)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
reg.form.SALL.plot
## `geom_smooth()` using formula = 'y ~ x'

reg.form.SBLL <- lm(form$SBLL ~ form$SL)
sd.form.SBLL <- rstandard(reg.form.SBLL)
reg.form.SBLL.plot <- ggplot(form, aes(x = SL, y = SBLL)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
reg.form.SBLL.plot
## `geom_smooth()` using formula = 'y ~ x'

reg.form.SBDF <- lm(form$SBDF ~ form$SL)
sd.form.SBDF <- rstandard(reg.form.SBDF)
reg.form.SBDF.plot <- ggplot(form, aes(x = SL, y = SBDF)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
reg.form.SBDF.plot
## `geom_smooth()` using formula = 'y ~ x'

reg.form.BD <- lm(form$BD_log ~ form$SL_log)
sd.form.BD <- rstandard(reg.form.BD)
reg.form.BD.plot <- ggplot(form, aes(x = SL_log, y = BD_log)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
reg.form.BD.plot
## `geom_smooth()` using formula = 'y ~ x'

reg.form.CPD <- lm(form$CPD_log ~ form$SL_log)
sd.form.CPD <- rstandard(reg.form.CPD)
reg.form.CPD.plot <- ggplot(form, aes(x = SL_log, y = CPD_log)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
reg.form.CPD.plot
## `geom_smooth()` using formula = 'y ~ x'

reg.form.CPL <- lm(form$CPL_log ~ form$SL_log)
sd.form.CPL <- rstandard(reg.form.CPL)
reg.form.CPL.plot <- ggplot(form, aes(x = SL_log, y = CPL_log)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
reg.form.CPL.plot
## `geom_smooth()` using formula = 'y ~ x'

reg.form.PreDL <- lm(form$PreDL_log ~ form$SL_log)
sd.form.PreDL <- rstandard(reg.form.PreDL)
reg.form.PreDL.plot <- ggplot(form, aes(x = SL_log, y = PreDL_log)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
reg.form.PreDL.plot
## `geom_smooth()` using formula = 'y ~ x'

reg.form.DbL <- lm(form$DbL_log ~ form$SL_log)
sd.form.DbL <- rstandard(reg.form.DbL)
reg.form.DbL.plot <- ggplot(form, aes(x = SL_log, y = DbL_log)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
reg.form.DbL.plot
## `geom_smooth()` using formula = 'y ~ x'

reg.form.HL <- lm(form$HL_log ~ form$SL_log)
sd.form.HL <- rstandard(reg.form.HL)
reg.form.HL.plot <- ggplot(form, aes(x = SL_log, y = HL_log)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
reg.form.HL.plot
## `geom_smooth()` using formula = 'y ~ x'

reg.form.HD <- lm(form$HD_log ~ form$SL_log)
sd.form.HD <- rstandard(reg.form.HD)
reg.form.HD.plot <- ggplot(form, aes(x = SL_log, y = HD_log)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
reg.form.HD.plot
## `geom_smooth()` using formula = 'y ~ x'

reg.form.HW <- lm(form$HW_log ~ form$SL_log)
sd.form.HW <- rstandard(reg.form.HW)
reg.form.HW.plot <- ggplot(form, aes(x = SL_log, y = HW_log)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
reg.form.HW.plot
## `geom_smooth()` using formula = 'y ~ x'

reg.form.SnL <- lm(form$SnL_log ~ form$SL_log)
sd.form.SnL <- rstandard(reg.form.SnL)
reg.form.SnL.plot <- ggplot(form, aes(x = SL_log, y = SnL_log)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
reg.form.SnL.plot
## `geom_smooth()` using formula = 'y ~ x'

reg.form.OL <- lm(form$OL_log ~ form$SL_log)
sd.form.OL <- rstandard(reg.form.OL)
reg.form.OL.plot <- ggplot(form, aes(x = SL_log, y = OL_log)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
reg.form.OL.plot
## `geom_smooth()` using formula = 'y ~ x'

reg.form.FLA <- lm(form$FLA ~ form$SL)
sd.form.FLA <- rstandard(reg.form.FLA)
reg.form.FLA.plot <- ggplot(form, aes(x = SL, y = FLA)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
reg.form.FLA.plot
## `geom_smooth()` using formula = 'y ~ x'

##### MEX #####
reg.mex.D <- lm(mex$D ~ mex$SL)
sd.mex.D <- rstandard(reg.mex.D)
reg.mex.D.plot <- ggplot(mex, aes(x = SL, y = D)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
reg.mex.D.plot
## `geom_smooth()` using formula = 'y ~ x'

reg.mex.P1 <- lm(mex$P1 ~ mex$SL)
sd.mex.P1 <- rstandard(reg.mex.P1)
reg.mex.P1.plot <- ggplot(mex, aes(x = SL, y = P1)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
reg.mex.P1.plot
## `geom_smooth()` using formula = 'y ~ x'

reg.mex.P2.L <- lm(mex$P2.L ~ mex$SL)
sd.mex.P2.L <- rstandard(reg.mex.P2.L)
reg.mex.P2.L.plot <- ggplot(mex, aes(x = SL, y = P2.L)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
reg.mex.P2.L.plot
## `geom_smooth()` using formula = 'y ~ x'

reg.mex.P2.R <- lm(mex$P2.R ~ mex$SL)
sd.mex.P2.R <- rstandard(reg.mex.P2.R)
reg.mex.P2.R.plot <- ggplot(mex, aes(x = SL, y = P2.R)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
reg.mex.P2.R.plot
## `geom_smooth()` using formula = 'y ~ x'

reg.mex.A <- lm(mex$A ~ mex$SL)
sd.mex.A <- rstandard(reg.mex.A)
reg.mex.A.plot <- ggplot(mex, aes(x = SL, y = A)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
reg.mex.A.plot
## `geom_smooth()` using formula = 'y ~ x'

reg.mex.P1.R <- lm(mex$P1.R ~ mex$SL)
sd.mex.P1.R <- rstandard(reg.mex.P1.R)
reg.mex.P1.R.plot <- ggplot(mex, aes(x = SL, y = P1.R)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
reg.mex.P1.R.plot
## `geom_smooth()` using formula = 'y ~ x'

reg.mex.LLSC <- lm(mex$LLSC ~ mex$SL)
sd.mex.LLSC <- rstandard(reg.mex.LLSC)
reg.mex.LLSC.plot <- ggplot(mex, aes(x = SL, y = LLSC)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
reg.mex.LLSC.plot
## `geom_smooth()` using formula = 'y ~ x'

reg.mex.SALL <- lm(mex$SALL ~ mex$SL)
sd.mex.SALL <- rstandard(reg.mex.SALL)
reg.mex.SALL.plot <- ggplot(mex, aes(x = SL, y = SALL)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
reg.mex.SALL.plot
## `geom_smooth()` using formula = 'y ~ x'

reg.mex.SBLL <- lm(mex$SBLL ~ mex$SL)
sd.mex.SBLL <- rstandard(reg.mex.SBLL)
reg.mex.SBLL.plot <- ggplot(mex, aes(x = SL, y = SBLL)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
reg.mex.SBLL.plot
## `geom_smooth()` using formula = 'y ~ x'

reg.mex.SBDF <- lm(mex$SBDF ~ mex$SL)
sd.mex.SBDF <- rstandard(reg.mex.SBDF)
reg.mex.SBDF.plot <- ggplot(mex, aes(x = SL, y = SBDF)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
reg.mex.SBDF.plot
## `geom_smooth()` using formula = 'y ~ x'

reg.mex.BD <- lm(mex$BD_log ~ mex$SL_log)
sd.mex.BD <- rstandard(reg.mex.BD)
reg.mex.BD.plot <- ggplot(mex, aes(x = SL_log, y = BD_log)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
reg.mex.BD.plot
## `geom_smooth()` using formula = 'y ~ x'

reg.mex.CPD <- lm(mex$CPD_log ~ mex$SL_log)
sd.mex.CPD <- rstandard(reg.mex.CPD)
reg.mex.CPD.plot <- ggplot(mex, aes(x = SL_log, y = CPD_log)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
reg.mex.CPD.plot
## `geom_smooth()` using formula = 'y ~ x'

reg.mex.CPL <- lm(mex$CPL_log ~ mex$SL_log)
sd.mex.CPL <- rstandard(reg.mex.CPL)
reg.mex.CPL.plot <- ggplot(mex, aes(x = SL_log, y = CPL_log)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
reg.mex.CPL.plot
## `geom_smooth()` using formula = 'y ~ x'

reg.mex.PreDL <- lm(mex$PreDL_log ~ mex$SL_log)
sd.mex.PreDL <- rstandard(reg.mex.PreDL)
reg.mex.PreDL.plot <- ggplot(mex, aes(x = SL_log, y = PreDL_log)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
reg.mex.PreDL.plot
## `geom_smooth()` using formula = 'y ~ x'

reg.mex.DbL <- lm(mex$DbL_log ~ mex$SL_log)
sd.mex.DbL <- rstandard(reg.mex.DbL)
reg.mex.DbL.plot <- ggplot(mex, aes(x = SL_log, y = DbL_log)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
reg.mex.DbL.plot
## `geom_smooth()` using formula = 'y ~ x'

reg.mex.HL <- lm(mex$HL_log ~ mex$SL_log)
sd.mex.HL <- rstandard(reg.mex.HL)
reg.mex.HL.plot <- ggplot(mex, aes(x = SL_log, y = HL_log)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
reg.mex.HL.plot
## `geom_smooth()` using formula = 'y ~ x'

reg.mex.HD <- lm(mex$HD_log ~ mex$SL_log)
sd.mex.HD <- rstandard(reg.mex.HD)
reg.mex.HD.plot <- ggplot(mex, aes(x = SL_log, y = HD_log)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
reg.mex.HD.plot
## `geom_smooth()` using formula = 'y ~ x'

reg.mex.HW <- lm(mex$HW_log ~ mex$SL_log)
sd.mex.HW <- rstandard(reg.mex.HW)
reg.mex.HW.plot <- ggplot(mex, aes(x = SL_log, y = HW_log)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
reg.mex.HW.plot
## `geom_smooth()` using formula = 'y ~ x'

reg.mex.SnL <- lm(mex$SnL_log ~ mex$SL_log)
sd.mex.SnL <- rstandard(reg.mex.SnL)
reg.mex.SnL.plot <- ggplot(mex, aes(x = SL_log, y = SnL_log)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
reg.mex.SnL.plot
## `geom_smooth()` using formula = 'y ~ x'

reg.mex.OL <- lm(mex$OL_log ~ mex$SL_log)
sd.mex.OL <- rstandard(reg.mex.OL)
reg.mex.OL.plot <- ggplot(mex, aes(x = SL_log, y = OL_log)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
reg.mex.OL.plot
## `geom_smooth()` using formula = 'y ~ x'

reg.mex.FLA <- lm(mex$FLA ~ mex$SL)
sd.mex.FLA <- rstandard(reg.mex.FLA)
reg.mex.FLA.plot <- ggplot(mex, aes(x = SL, y = FLA)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
reg.mex.FLA.plot
## `geom_smooth()` using formula = 'y ~ x'

Residuals

  • STEP TWO: get residuals for each individual for traits that were influenced by body size

  • STEP THREE: convert residuals to absolute value

##### LAT #####

abs.lat.D <- abs(res.lat.D)
mean(abs.lat.D)
## [1] 0.5375505
abs.lat.P1 <- abs(res.lat.P1)
mean(abs.lat.P1)
## [1] 0.5466667
abs.lat.P1.R <- abs(res.lat.P1.R)
mean(abs.lat.P1.R)
## [1] 0.5799393
abs.lat.LLSC <- abs(res.lat.LLSC)
mean(abs.lat.LLSC)
## [1] 0.7664044
abs.lat.SBLL <- abs(res.lat.SBLL)
mean(abs.lat.SBLL)
## [1] 0.2461336
abs.lat.BD <- abs(res.lat.BD)
mean(abs.lat.BD)
## [1] 0.03636357
abs.lat.CPD <- abs(res.lat.CPD)
mean(abs.lat.CPD)
## [1] 0.03695739
abs.lat.CPL <- abs(res.lat.CPL)
mean(abs.lat.CPL)
## [1] 0.03732998
abs.lat.PreDL <- abs(res.lat.PreDL)
mean(abs.lat.PreDL)
## [1] 0.04810736
abs.lat.DbL <- abs(res.lat.DbL)
mean(abs.lat.DbL)
## [1] 0.04126962
abs.lat.HL <- abs(res.lat.HL)
mean(abs.lat.HL)
## [1] 0.03588147
abs.lat.HD <- abs(res.lat.HD)
mean(abs.lat.HD)
## [1] 0.03506907
abs.lat.HW <- abs(res.lat.HW)
mean(abs.lat.HW)
## [1] 0.03208151
abs.lat.SnL <- abs(res.lat.SnL)
mean(abs.lat.SnL)
## [1] 0.03208694
abs.lat.OL <- abs(res.lat.OL)
mean(abs.lat.OL)
## [1] 1.416081e-17
##### FORM #####

abs.form.D <- abs(res.form.D)
mean(abs.form.D)
## [1] 0.5668177
abs.form.P1 <- abs(res.form.P1)
mean(abs.form.P1)
## [1] 0.4843616
abs.form.P1.R <- abs(res.form.P1.R)
mean(abs.form.P1.R)
## [1] 0.4242033
abs.form.LLSC <- abs(res.form.LLSC)
mean(abs.form.LLSC)
## [1] 0.9038801
abs.form.SBLL <- abs(res.form.SBLL)
mean(abs.form.SBLL)
## [1] 0.3399272
abs.form.BD <- abs(res.form.BD)
mean(abs.form.BD)
## [1] 0.04191086
abs.form.CPD <- abs(res.form.CPD)
mean(abs.form.CPD)
## [1] 0.03046568
abs.form.CPL <- abs(res.form.CPL)
mean(abs.form.CPL)
## [1] 0.03738523
abs.form.PreDL <- abs(res.form.PreDL)
mean(abs.form.PreDL)
## [1] 0.05207252
abs.form.DbL <- abs(res.form.DbL)
mean(abs.form.DbL)
## [1] 0.03261267
abs.form.HL <- abs(res.form.HL)
mean(abs.form.HL)
## [1] 0.03919057
abs.form.HD <- abs(res.form.HD)
mean(abs.form.HD)
## [1] 0.03507274
abs.form.HW <- abs(res.form.HW)
mean(abs.form.HW)
## [1] 0.03744209
abs.form.SnL <- abs(res.form.SnL)
mean(abs.form.SnL)
## [1] 0.03321741
abs.form.OL <- abs(res.form.OL)
mean(abs.form.OL)
## [1] 6.018764e-18
##### MEX #####

abs.mex.D <- abs(res.mex.D)
mean(abs.mex.D)
## [1] 0.1657275
abs.mex.P1 <- abs(res.mex.P1)
mean(abs.mex.P1)
## [1] 0.5686425
abs.mex.P1.R <- abs(res.mex.P1.R)
mean(abs.mex.P1.R)
## [1] 0.458723
abs.mex.LLSC <- abs(res.mex.LLSC)
mean(abs.mex.LLSC)
## [1] 0.4434954
abs.mex.SBLL <- abs(res.mex.SBLL)
mean(abs.mex.SBLL)
## [1] 0.2713231
abs.mex.BD <- abs(res.mex.BD)
mean(abs.mex.BD)
## [1] 0.03300381
abs.mex.CPD <- abs(res.mex.CPD)
mean(abs.mex.CPD)
## [1] 0.03277006
abs.mex.CPL <- abs(res.mex.CPL)
mean(abs.mex.CPL)
## [1] 0.03878753
abs.mex.PreDL <- abs(res.mex.PreDL)
mean(abs.mex.PreDL)
## [1] 0.03943373
abs.mex.DbL <- abs(res.mex.DbL)
mean(abs.mex.DbL)
## [1] 0.02950105
abs.mex.HL <- abs(res.mex.HL)
mean(abs.mex.HL)
## [1] 0.05562721
abs.mex.HD <- abs(res.mex.HD)
mean(abs.mex.HD)
## [1] 0.03668815
abs.mex.HW <- abs(res.mex.HW)
mean(abs.mex.HW)
## [1] 0.03882769
abs.mex.SnL <- abs(res.mex.SnL)
mean(abs.mex.SnL)
## [1] 0.04900589
abs.mex.OL <- abs(res.mex.OL)
mean(abs.mex.OL)
## [1] 4.001035e-18
#let's get this into the raw1 data set so that we can plot this more easily

abs.res.D <- c(abs.lat.D, abs.form.D, abs.mex.D)
abs.res.P1 <- c(abs.lat.P1, abs.form.P1, abs.mex.P1)
abs.res.P1.R <- c(abs.lat.P1.R, abs.form.P1.R, abs.mex.P1.R)
abs.res.LLSC<- c(abs.lat.LLSC, abs.form.LLSC, abs.mex.LLSC)
abs.res.SBLL<- c(abs.lat.SBLL, abs.form.SBLL, abs.mex.SBLL)
abs.res.BD<- c(abs.lat.BD, abs.form.BD, abs.mex.BD)
abs.res.CPD<- c(abs.lat.CPD, abs.form.CPD, abs.mex.CPD)
abs.res.CPL<- c(abs.lat.CPL, abs.form.CPL, abs.mex.CPL)
abs.res.PreDL <- c(abs.lat.PreDL, abs.form.PreDL, abs.mex.PreDL)
abs.res.DbL <- c(abs.lat.DbL, abs.form.DbL, abs.mex.DbL)
abs.res.HL<- c(abs.lat.HL, abs.form.HL, abs.mex.HL)
abs.res.HD<- c(abs.lat.HD, abs.form.HD, abs.mex.HD)
abs.res.HW <- c(abs.lat.HW, abs.form.HW, abs.mex.HW)
abs.res.SnL <- c(abs.lat.SnL, abs.form.SnL, abs.mex.SnL)
abs.res.OL <- c(abs.lat.OL, abs.form.OL, abs.mex.OL)


raw2 <- cbind(raw1, abs.res.D, abs.res.P1, abs.res.P1.R, abs.res.LLSC, abs.res.SBLL, abs.res.BD, abs.res.CPD, abs.res.CPL, abs.res.PreDL, abs.res.DbL, abs.res.HL, abs.res.HD, abs.res.HW, abs.res.SnL, abs.res.OL)

Residual Comparisons

  • STEP FOUR: plot the mean of the |residuals| for both species, for all traits influenced by body size
library(ggbeeswarm)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
ggplot(raw2, aes(SPP, abs.res.D)) +
  geom_point(alpha=0.3) +
  stat_summary(fun.data=function(x){mean_cl_normal(x, conf.int=.683)}, geom="errorbar", 
               width=0.03, colour="red", alpha=0.7) +
  stat_summary(fun=mean, geom="point", fill="red", pch=21, size=3)

scatter_violin <- ggplot(data=raw2, aes(x=SPP, y=abs.res.D)) +
  geom_violin(trim = FALSE) + 
  stat_summary(
    fun.data = "mean_sdl",  fun.args = list(mult = 1),
    geom = "pointrange", color = "black"
    )

print(scatter_violin)

scatter_violin1 <- ggplot(data=raw2, aes(x=SPP, y=abs.res.D)) +
  geom_violin(trim = FALSE) + 
  stat_summary(fun.data=function(x){mean_cl_normal(x, conf.int=.683)}, geom="errorbar", 
               width=0.03, colour="red", alpha=0.7) +
  stat_summary(fun=mean, geom="crossbar", fill="red", width=0.03)

print(scatter_violin1)

ggplot(raw2, aes(SPP, abs.res.P1)) +
  geom_point(alpha=0.3) +
  stat_summary(fun.data=function(x){mean_cl_normal(x, conf.int=.683)}, geom="errorbar", 
               width=0.03, colour="red", alpha=0.7) +
  stat_summary(fun=mean, geom="point", fill="red", pch=21, size=3)

ggplot(raw2, aes(SPP, abs.res.P1.R)) +
  geom_point(alpha=0.3) +
  stat_summary(fun.data=function(x){mean_cl_normal(x, conf.int=.683)}, geom="errorbar", 
               width=0.03, colour="red", alpha=0.7) +
  stat_summary(fun=mean, geom="point", fill="red", pch=21, size=3) 

ggplot(raw2, aes(SPP, abs.res.LLSC)) +
  geom_point(alpha=0.3) +
  stat_summary(fun.data=function(x){mean_cl_normal(x, conf.int=.683)}, geom="errorbar", 
               width=0.03, colour="red", alpha=0.7) +
  stat_summary(fun=mean, geom="point", fill="red", pch=21, size=3)

ggplot(raw2, aes(SPP, abs.res.SBLL)) +
  geom_point(alpha=0.3) +
  stat_summary(fun.data=function(x){mean_cl_normal(x, conf.int=.683)}, geom="errorbar", 
               width=0.03, colour="red", alpha=0.7) +
  stat_summary(fun=mean, geom="point", fill="red", pch=21, size=3)

ggplot(raw2, aes(SPP, abs.res.BD)) +
  geom_point(alpha=0.3) +
  stat_summary(fun.data=function(x){mean_cl_normal(x, conf.int=.683)}, geom="errorbar", 
               width=0.03, colour="red", alpha=0.7) +
  stat_summary(fun=mean, geom="point", fill="red", pch=21, size=3)

ggplot(raw2, aes(SPP, abs.res.CPD)) +
  geom_point(alpha=0.3) +
  stat_summary(fun.data=function(x){mean_cl_normal(x, conf.int=.683)}, geom="errorbar", 
               width=0.03, colour="red", alpha=0.7) +
  stat_summary(fun=mean, geom="point", fill="red", pch=21, size=3)

ggplot(raw2, aes(SPP, abs.res.CPL)) +
  geom_point(alpha=0.3) +
  stat_summary(fun.data=function(x){mean_cl_normal(x, conf.int=.683)}, geom="errorbar", 
               width=0.03, colour="red", alpha=0.7) +
  stat_summary(fun=mean, geom="point", fill="red", pch=21, size=3) 

ggplot(raw2, aes(SPP, abs.res.PreDL)) +
  geom_point(alpha=0.3) +
  stat_summary(fun.data=function(x){mean_cl_normal(x, conf.int=.683)}, geom="errorbar", 
               width=0.03, colour="red", alpha=0.7) +
  stat_summary(fun=mean, geom="point", fill="red", pch=21, size=3)

ggplot(raw2, aes(SPP, abs.res.DbL)) +
  geom_point(alpha=0.3) +
  stat_summary(fun.data=function(x){mean_cl_normal(x, conf.int=.683)}, geom="errorbar", 
               width=0.03, colour="red", alpha=0.7) +
  stat_summary(fun=mean, geom="point", fill="red", pch=21, size=3)

ggplot(raw2, aes(SPP, abs.res.HL)) +
  geom_point(alpha=0.3) +
  stat_summary(fun.data=function(x){mean_cl_normal(x, conf.int=.683)}, geom="errorbar", 
               width=0.03, colour="red", alpha=0.7) +
  stat_summary(fun=mean, geom="point", fill="red", pch=21, size=3)

ggplot(raw2, aes(SPP, abs.res.HD)) +
  geom_point(alpha=0.3) +
  stat_summary(fun.data=function(x){mean_cl_normal(x, conf.int=.683)}, geom="errorbar", 
               width=0.03, colour="red", alpha=0.7) +
  stat_summary(fun=mean, geom="point", fill="red", pch=21, size=3)

ggplot(raw2, aes(SPP, abs.res.HW)) +
  geom_point(alpha=0.3) +
  stat_summary(fun.data=function(x){mean_cl_normal(x, conf.int=.683)}, geom="errorbar", 
               width=0.03, colour="red", alpha=0.7) +
  stat_summary(fun=mean, geom="point", fill="red", pch=21, size=3)

ggplot(raw2, aes(SPP, abs.res.SnL)) +
  geom_point(alpha=0.3) +
  stat_summary(fun.data=function(x){mean_cl_normal(x, conf.int=.683)}, geom="errorbar", 
               width=0.03, colour="red", alpha=0.7) +
  stat_summary(fun=mean, geom="point", fill="red", pch=21, size=3)

ggplot(raw2, aes(SPP, abs.res.OL)) +
  geom_point(alpha=0.3) +
  stat_summary(fun.data=function(x){mean_cl_normal(x, conf.int=.683)}, geom="errorbar", 
               width=0.03, colour="red", alpha=0.7) +
  stat_summary(fun=mean, geom="point", fill="red", pch=21, size=3)

Comparing variation

Variance tests

t-test for residuals, continuous f-test for raw, continuous all discrete will go through glmer/non parametric with either residuals or raw

(note: did non-parametric instead of transforming, but Mike said to transform, so I copied that work to the ‘morphology-analysis-final’ Rmd).

F-test on residuals doesn’t make much sense; the residuals themselves are representative of the variation present, as they are the distance from the mean. Therefore, F-test on residuals is like variance of the variance. Instead, I have to do a t-test on the absolute value of the residuals. In this sense, we want to see if the mean of the absolute residuals is higher or lower for asexual species–is the average amount of variation higher or lower for this trait? Based on the regressions, if the trait was influenced by body size, I will perform a t-test on the absolute value of the residuals. If the trait was not influenced by body size, I will perform an F-test of variance on the raw data.

Quick results summary: For the F-test on raw data, none of the traits were significantly different (P2L/R, A, SBDF, FLA). For the t-tests on residuals, the only significant traits are left pectoral (), right pectoral (lat>form), scales above lateral line (), scales below lateral line (form>lat), and head length ().

-   will do two-tailed and check out the residual means to infer direction; for traits in which we use raw data, a one-tailed f-test will be perfomed in both direction to determine which species is varying more. We will also visulize the variation using a histogram to confirm direction results.

F tests to compare species

ALL OF THESE ARE DISCRETE SO I DON’T THINK I DO ANY F TESTS

This will only be for the traits that did not have a significant regression compared to SL (so P2s, A, SALL, SBDF, and FLA); none of these traits were transformed, so just doing raw for now.

var.left.pel <- var.test(lat\(P2.L, form\)P2.L, alternative = “two.sided”) var.left.pel

var.rig.pel <- var.test(lat\(P2.R, form\)P2.R, alternative = “two.sided”) var.rig.pel

var.anal <- var.test(lat\(A, form\)A, alternative = “two.sided”) var.anal

var.sca.ab.ll <- var.test(lat\(SALL, form\)SALL, alternative = “two.sided”) var.sca.ab.ll

var.sca.df <- var.test(lat\(SBDF, form\)SBDF, alternative = “two.sided”) var.sca.df

var.flu.asy <- var.test(lat\(FLA, form\)FLA, alternative = “two.sided”) var.flu.asy





T tests to compare species

CONTINUOUS RESIDUALS

  • will do two-tailed and check out the means to infer direction (see beginning for trait means w/o body correction)
  • this will be for the traits that did have a significant regression compared to SL; will be using the absolute value residuals.
  • For continuous traits, we will be using the log transformed data.
## 
##  Welch Two Sample t-test
## 
## data:  abs.lat.BD and abs.form.BD
## t = -0.83822, df = 215.6, p-value = 0.4028
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.018591485  0.007496913
## sample estimates:
##  mean of x  mean of y 
## 0.03636357 0.04191086
## 
##  Welch Two Sample t-test
## 
## data:  abs.lat.CPD and abs.form.CPD
## t = 2.1468, df = 282.5, p-value = 0.03266
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  0.0005394989 0.0124439196
## sample estimates:
##  mean of x  mean of y 
## 0.03695739 0.03046568
## 
##  Welch Two Sample t-test
## 
## data:  abs.lat.CPL and abs.form.CPL
## t = -0.015652, df = 291.31, p-value = 0.9875
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.007002590  0.006892088
## sample estimates:
##  mean of x  mean of y 
## 0.03732998 0.03738523
## 
##  Welch Two Sample t-test
## 
## data:  abs.lat.PreDL and abs.form.PreDL
## t = -0.67807, df = 215.55, p-value = 0.4985
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.015491199  0.007560893
## sample estimates:
##  mean of x  mean of y 
## 0.04810736 0.05207252
## 
##  Welch Two Sample t-test
## 
## data:  abs.lat.DbL and abs.form.DbL
## t = 2.3458, df = 281.07, p-value = 0.01968
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  0.00139255 0.01592134
## sample estimates:
##  mean of x  mean of y 
## 0.04126962 0.03261267
## 
##  Welch Two Sample t-test
## 
## data:  abs.lat.HL and abs.form.HL
## t = -0.90175, df = 296.34, p-value = 0.3679
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.010530987  0.003912795
## sample estimates:
##  mean of x  mean of y 
## 0.03588147 0.03919057
## 
##  Welch Two Sample t-test
## 
## data:  abs.lat.HD and abs.form.HD
## t = -0.00099329, df = 267.38, p-value = 0.9992
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.007275651  0.007268313
## sample estimates:
##  mean of x  mean of y 
## 0.03506907 0.03507274
## 
##  Welch Two Sample t-test
## 
## data:  abs.lat.HW and abs.form.HW
## t = -1.6871, df = 297.98, p-value = 0.09263
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.0116134871  0.0008923157
## sample estimates:
##  mean of x  mean of y 
## 0.03208151 0.03744209
## 
##  Welch Two Sample t-test
## 
## data:  abs.lat.SnL and abs.form.SnL
## t = -0.38886, df = 295.99, p-value = 0.6977
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -0.006851803  0.004590862
## sample estimates:
##  mean of x  mean of y 
## 0.03208694 0.03321741
## 
##  Welch Two Sample t-test
## 
## data:  abs.lat.OL and abs.form.OL
## t = 1.318, df = 152.51, p-value = 0.1895
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -4.062607e-18  2.034670e-17
## sample estimates:
##    mean of x    mean of y 
## 1.416081e-17 6.018764e-18

Levene’s test

Since the count data can’t do parametric, I will do non-parametric tests. Logan mentioned a glmer with poisson distribution but I don’t think this will answer the question of a difference in variance between the species. The linear model will instead tell me how much species effects differences in the trait?? I guess?

I will use Levene’s test as a non-parametric F-test, to check out differences in variance in raw data.

I will use Mann Whitney U test as a non-parametric t-test, to check out differences in variance in residual data.

Still only performing this on the traits that did NOT vary with SL (P2L/R, A, SALL, SBLL [between sail and amz only], SBDF, FLA).

raw3 <- raw2[raw2$SPP !="p.mexicana", ]


library(car)
## Loading required package: carData
## 
## Attaching package: 'car'
## The following object is masked from 'package:dplyr':
## 
##     recode
library(carData)

(LT_P2L <- leveneTest(P2.L~SPP, data=raw3)) #gives nothing since it's all the same value
## Levene's Test for Homogeneity of Variance (center = median)
##        Df F value Pr(>F)
## group   1     NaN    NaN
##       298
(LT_P2R <- leveneTest(P2.R~SPP, data=raw3)) #same as above
## Levene's Test for Homogeneity of Variance (center = median)
##        Df F value Pr(>F)
## group   1     NaN    NaN
##       298
(LT_A <- leveneTest(A~SPP, data=raw3)) 
## Levene's Test for Homogeneity of Variance (center = median)
##        Df F value Pr(>F)
## group   1  0.3962 0.5295
##       298
(LT_SALL <- leveneTest(SALL~SPP, data=raw3))
## Levene's Test for Homogeneity of Variance (center = median)
##        Df F value Pr(>F)
## group   1  0.8627 0.3537
##       298
(LT_SBLL <- leveneTest(SBLL~SPP, data=raw3))
## Levene's Test for Homogeneity of Variance (center = median)
##        Df F value Pr(>F)  
## group   1    3.36 0.0678 .
##       298                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(LT_SBDF <- leveneTest(SBDF~SPP, data=raw3))
## Levene's Test for Homogeneity of Variance (center = median)
##        Df F value Pr(>F)
## group   1  0.0081 0.9283
##       298
(LT_FLA <- leveneTest(FLA~SPP, data=raw3))
## Levene's Test for Homogeneity of Variance (center = median)
##        Df F value Pr(>F)
## group   1  2.7164 0.1004
##       298

Mann Whitney U tests

This will be performed on traits that DID vary with SL.

(MW_D <- wilcox.test(abs.res.D~SPP, data=raw3, conf.int=T))
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  abs.res.D by SPP
## W = 12186, p-value = 0.1545
## alternative hypothesis: true location shift is not equal to 0
## 95 percent confidence interval:
##  -0.05689062  0.19840357
## sample estimates:
## difference in location 
##              0.1030553
(MW_P1 <- wilcox.test(abs.res.P1~SPP, data=raw3, conf.int=T))
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  abs.res.P1 by SPP
## W = 8765, p-value = 0.001606
## alternative hypothesis: true location shift is not equal to 0
## 95 percent confidence interval:
##  -0.23328614 -0.07245187
## sample estimates:
## difference in location 
##             -0.1618195
(MW_P1R <- wilcox.test(abs.res.P1.R~SPP, data=raw3, conf.int=T))
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  abs.res.P1.R by SPP
## W = 7363, p-value = 4.86e-07
## alternative hypothesis: true location shift is not equal to 0
## 95 percent confidence interval:
##  -0.2596445 -0.1543936
## sample estimates:
## difference in location 
##             -0.2082393
(MW_LLSC <- wilcox.test(abs.res.LLSC~SPP, data=raw3, conf.int=T))
## 
##  Wilcoxon rank sum test with continuity correction
## 
## data:  abs.res.LLSC by SPP
## W = 12119, p-value = 0.1822
## alternative hypothesis: true location shift is not equal to 0
## 95 percent confidence interval:
##  -0.03418245  0.19226705
## sample estimates:
## difference in location 
##              0.0763304

Summary of variance results

average.table <- matrix(c(MW_D\(p.value, mean(abs.lat.D, na.rm = TRUE), mean(abs.form.D, na.rm = TRUE), MW_P1\)p.value, mean(abs.lat.P1, na.rm = TRUE), mean(abs.form.P1, na.rm = TRUE), LT_P2R\("Pr(>F)", mean(lat\)P2.R, na.rm = TRUE), mean(form\(P2.R, na.rm = TRUE), LT_P2R\)“Pr(>F)”, mean(lat\(P2.L, na.rm = TRUE), mean(form\)P2.L, na.rm = TRUE), LT_A\("Pr(>F)", mean(lat\)A, na.rm = TRUE), mean(form\(A, na.rm = TRUE), MW_P1R\)p.value, mean(abs.lat.P1.R, na.rm = TRUE), mean(abs.form.P1.R, na.rm = TRUE), MW_LLSC\(p.value, mean(abs.lat.LLSC, na.rm = TRUE), mean(abs.form.LLSC, na.rm = TRUE), LT_SALL\)“Pr(>F)”, mean(lat\(SALL, na.rm = TRUE), mean(form\)SALL, na.rm = TRUE), LT_SBLL\("Pr(>F)", mean(abs.lat.SBLL, na.rm = TRUE), mean(abs.form.SBLL, na.rm = TRUE), LT_SBDF\)“Pr(>F)”, mean(lat\(SBDF, na.rm = TRUE), mean(form\)SBDF, na.rm = TRUE), ttest.abs.BD\(p.value, mean(abs.lat.BD, na.rm = TRUE), mean(abs.form.BD, na.rm = TRUE), ttest.abs.CPD\)p.value, mean(abs.lat.CPD, na.rm = TRUE), mean(abs.form.CPD, na.rm = TRUE), ttest.abs.CPL\(p.value, mean(abs.lat.CPL, na.rm = TRUE), mean(abs.form.CPL, na.rm = TRUE), ttest.abs.PreDL\)p.value, mean(abs.lat.PreDL, na.rm = TRUE), mean(abs.form.PreDL, na.rm = TRUE), ttest.abs.DbL\(p.value, mean(abs.lat.DbL, na.rm = TRUE), mean(abs.form.DbL, na.rm = TRUE), ttest.abs.HL\)p.value, mean(abs.lat.HL, na.rm = TRUE), mean(abs.form.HL, na.rm = TRUE), ttest.abs.HD\(p.value, mean(abs.lat.HD, na.rm = TRUE), mean(abs.form.HD, na.rm = TRUE), ttest.abs.HW\)p.value, mean(abs.lat.HW, na.rm = TRUE), mean(abs.form.HW, na.rm = TRUE), ttest.abs.SnL\(p.value, mean(abs.lat.SnL, na.rm = TRUE), mean(abs.form.SnL, na.rm = TRUE), ttest.abs.OL\)p.value, mean(abs.lat.OL, na.rm = TRUE), mean(abs.form.OL, na.rm = TRUE), LT_FLA\("Pr(>F)", mean(lat\)FLA, na.rm = TRUE), mean(form$FLA, na.rm = TRUE)), ncol=3, byrow=TRUE)

colnames(average.table)<- c(“p-value for variance”, “lat mean”, “form mean”)

rownames(average.table) <- c(“dorsal rays”, “L pect rays”, “R pelvic rays”, ”L pelvic rays”, “Anal rays”, “R pect rays”, ”lat line scales”, ”scales above LL”, ”scales below LL”, ”scales before dor”, “body dep”, “peduncle dep”, “peduncle leng”, “pre-dorsal”, “dorsal base”, “head length”, “head depth”, “head width”, “snout leng”, “orbital leng”, “fluct. asymmetries”)

average.table

Variance plots

plot_multi_histogram <- function(df, feature, label_column) {
    plt <- ggplot(df, aes(x=eval(parse(text=feature)), fill=eval(parse(text=label_column)))) +
    geom_histogram(bins=10, alpha=0.4, position="identity", aes(y = ..density..), color="black") +
    geom_density(alpha=0.2)  +
    labs(x=feature, y = "Density")
    plt + guides(fill=guide_legend(title=label_column))
}

plot_multi_histogram(raw3, "abs.res.CPD", "SPP")
## Warning: The dot-dot notation (`..density..`) was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(density)` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

plot_multi_histogram(raw3, "abs.res.DbL", "SPP")

plot_multi_histogram(raw3, "abs.res.P1", "SPP")

plot_multi_histogram(raw3, "abs.res.P1.R", "SPP")

ANOVA

Will run an ANOVA on the residuals with location and species as fixed effects. This will show me if morphology depends on the species, the location, and if the location and species interact to determine morphology. CANNOT DO WITH COUNT DATA–may have to run a glmer with zero-tuncated negative binomial/quasi-poisson distribution instead (think this will tell me if species, location, and interaction affect the trait…). Will only work with the residuals, because comparing the means of the residuals indicates differences in variation. Would I not just use Levene’s and MWU again but include spp and zone/basin/watershed?

Leaving the anovas for count data in rn until I figure out what to do.

I will first run this using the zones as the location factor. Zones (1-4) represent the latitude range with equivalent sample sizes in each, since the collections were not equally representative of all latitudes, and I wanted to avoid a sampling bias when randomly selecting samples. Zone 1 corresponds to the southern most latitude range, and zone 4 corresponds to the northern most latitude range.

I will then run the same analysis using basin as the location factor. Since fish are physically isolated to the river basins they occupy, the genetic variation is also limited to that basin. Thus it is possible for fish within the same basin to be more similar due to genetic/physical constraints. (will also do with watershed just to see).

Lastly I will run ANOVAs with both zones and basins but with standardized residuals. This would allow me to compare overall variation across traits (at least those that are depended on body size) rather than just one trait at a time. Not 100% sure if this is useful (or correct to do), but thought it would be interesting.

ANOVAs

Amazon

Zones
library(ggplot2)

lat3 <- raw2[raw2$SPP == "p.latipinna", ]

form3 <- raw2[raw2$SPP == "p.latipinna", ]

A.D <- aov(abs.res.D ~ QUARTILE, data=form3)
summary(A.D)
##              Df Sum Sq Mean Sq F value  Pr(>F)   
## QUARTILE      3  2.306  0.7685   4.292 0.00635 **
## Residuals   130 23.277  0.1791                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ggplot(form3, aes(x=factor(QUARTILE), y=abs.res.D)) +
  geom_boxplot()

A.P1 <- aov(abs.res.P1 ~ QUARTILE, data=form3)
summary(A.P1)
##              Df Sum Sq Mean Sq F value Pr(>F)
## QUARTILE      3  0.159 0.05296   0.444  0.722
## Residuals   130 15.511 0.11932
A.P1.R <- aov(abs.res.P1.R ~ QUARTILE, data=form3)
summary(A.P1.R)
##              Df Sum Sq Mean Sq F value Pr(>F)
## QUARTILE      3  0.065 0.02172   0.166  0.919
## Residuals   130 16.987 0.13067
ggplot(form3, aes(x=factor(QUARTILE), y=abs.res.P1.R)) +
  geom_boxplot()

A.LLSC <- aov(abs.res.LLSC ~ QUARTILE, data=form3)
summary(A.LLSC)
##              Df Sum Sq Mean Sq F value Pr(>F)
## QUARTILE      3   1.31  0.4354   1.209  0.309
## Residuals   130  46.80  0.3600
A.SBLL <- aov(abs.res.SBLL ~ QUARTILE, data=form3)
summary(A.SBLL)
##              Df Sum Sq Mean Sq F value Pr(>F)
## QUARTILE      3  0.436  0.1454   1.142  0.335
## Residuals   130 16.546  0.1273
ggplot(form3, aes(x=factor(QUARTILE), y=abs.res.SBLL)) +
  geom_boxplot()

A.BD <- aov(abs.res.BD ~ QUARTILE, data=form3)
summary(A.BD)
##              Df  Sum Sq   Mean Sq F value Pr(>F)
## QUARTILE      3 0.00259 0.0008622   1.058  0.369
## Residuals   130 0.10591 0.0008147
ggplot(form3, aes(x=factor(QUARTILE), y=abs.res.BD)) +
  geom_boxplot()

A.CPD <- aov(abs.res.CPD ~ QUARTILE, data=form3)
summary(A.CPD)
##              Df  Sum Sq  Mean Sq F value  Pr(>F)    
## QUARTILE      3 0.01249 0.004163   6.832 0.00026 ***
## Residuals   130 0.07922 0.000609                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ggplot(form3, aes(x=factor(QUARTILE), y=abs.res.CPD)) +
  geom_boxplot()

A.CPL <- aov(abs.res.CPL ~ QUARTILE, data=form3)
summary(A.CPL)
##              Df  Sum Sq  Mean Sq F value   Pr(>F)    
## QUARTILE      3 0.02074 0.006912   9.448 1.08e-05 ***
## Residuals   130 0.09511 0.000732                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ggplot(form3, aes(x=factor(QUARTILE), y=abs.res.CPL)) +
  geom_boxplot()

A.PreDL <- aov(abs.res.PreDL ~ QUARTILE, data=form3)
summary(A.PreDL)
##              Df Sum Sq  Mean Sq F value Pr(>F)
## QUARTILE      3 0.0082 0.002742    0.79  0.502
## Residuals   130 0.4511 0.003470
ggplot(form3, aes(x=factor(QUARTILE), y=abs.res.PreDL)) +
  geom_boxplot()

A.DbL <- aov(abs.res.DbL ~ QUARTILE, data=form3)
summary(A.DbL)
##              Df Sum Sq  Mean Sq F value   Pr(>F)    
## QUARTILE      3 0.0195 0.006501   7.136 0.000179 ***
## Residuals   130 0.1184 0.000911                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ggplot(form3, aes(x=factor(QUARTILE), y=abs.res.DbL)) +
  geom_boxplot()

A.HL <- aov(abs.res.HL ~ QUARTILE, data=form3)
summary(A.HL)
##              Df  Sum Sq   Mean Sq F value Pr(>F)
## QUARTILE      3 0.00344 0.0011466   1.324  0.269
## Residuals   130 0.11260 0.0008661
ggplot(form3, aes(x=factor(QUARTILE), y=abs.res.HL)) +
  geom_boxplot()

A.HD <- aov(abs.res.HD ~ QUARTILE, data=form3)
summary(A.HD)
##              Df  Sum Sq  Mean Sq F value Pr(>F)
## QUARTILE      3 0.00672 0.002239    2.04  0.111
## Residuals   130 0.14269 0.001098
A.HW <- aov(abs.res.HW ~ QUARTILE, data=form3)
summary(A.HW)
##              Df  Sum Sq  Mean Sq F value   Pr(>F)    
## QUARTILE      3 0.01136 0.003788   7.212 0.000163 ***
## Residuals   130 0.06827 0.000525                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ggplot(form3, aes(x=factor(QUARTILE), y=abs.res.HW)) +
  geom_boxplot()

A.SnL <- aov(abs.res.SnL ~ QUARTILE, data=form3)
summary(A.SnL)
##              Df  Sum Sq   Mean Sq F value  Pr(>F)   
## QUARTILE      3 0.00766 0.0025545   5.052 0.00242 **
## Residuals   130 0.06573 0.0005056                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
A.OL <- aov(abs.res.OL ~ QUARTILE, data=form3)
summary(A.OL)
##              Df   Sum Sq   Mean Sq F value Pr(>F)
## QUARTILE      3 2.17e-32 7.242e-33   1.538  0.208
## Residuals   130 6.12e-31 4.708e-33
ggplot(form3, aes(x=factor(QUARTILE), y=abs.res.OL)) +
  geom_boxplot()

Basins
A1.D <- aov(abs.res.D ~ BASIN, data=form3)
summary(A1.D)
##              Df Sum Sq Mean Sq F value Pr(>F)  
## BASIN         5  2.789  0.5579   3.133 0.0106 *
## Residuals   128 22.793  0.1781                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ggplot(form3, aes(x=factor(BASIN), y=abs.res.D)) +
  geom_boxplot()

A1.P1 <- aov(abs.res.P1 ~ BASIN, data=form3)
summary(A1.P1)
##              Df Sum Sq Mean Sq F value Pr(>F)
## BASIN         5  0.945  0.1891   1.643  0.153
## Residuals   128 14.725  0.1150
A1.P1.R <- aov(abs.res.P1.R ~ BASIN, data=form3)
summary(A1.P1.R)
##              Df Sum Sq Mean Sq F value Pr(>F)
## BASIN         5  0.159 0.03187   0.241  0.943
## Residuals   128 16.893 0.13197
ggplot(form3, aes(x=factor(BASIN), y=abs.res.P1.R)) +
  geom_boxplot()

A1.LLSC <- aov(abs.res.LLSC ~ BASIN, data=form3)
summary(A1.LLSC)
##              Df Sum Sq Mean Sq F value Pr(>F)  
## BASIN         5   4.54  0.9080   2.668  0.025 *
## Residuals   128  43.57  0.3404                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ggplot(form3, aes(x=factor(BASIN), y=abs.res.LLSC)) +
  geom_boxplot()

A1.SBLL <- aov(abs.res.SBLL ~ BASIN, data=form3)
summary(A1.SBLL)
##              Df Sum Sq Mean Sq F value Pr(>F)
## BASIN         5  0.755  0.1510   1.191  0.317
## Residuals   128 16.228  0.1268
ggplot(form3, aes(x=factor(BASIN), y=abs.res.SBLL)) +
  geom_boxplot()

A1.BD <- aov(abs.res.BD ~ BASIN, data=form3)
summary(A1.BD)
##              Df  Sum Sq   Mean Sq F value Pr(>F)
## BASIN         5 0.00316 0.0006313   0.767  0.575
## Residuals   128 0.10534 0.0008230
ggplot(form3, aes(x=factor(BASIN), y=abs.res.BD)) +
  geom_boxplot()

A1.CPD <- aov(abs.res.CPD ~ BASIN, data=form3)
summary(A1.CPD)
##              Df  Sum Sq  Mean Sq F value   Pr(>F)    
## BASIN         5 0.01713 0.003425   5.879 6.34e-05 ***
## Residuals   128 0.07458 0.000583                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ggplot(form3, aes(x=factor(BASIN), y=abs.res.CPD)) +
  geom_boxplot()

A1.CPL <- aov(abs.res.CPL ~ BASIN, data=form3)
summary(A1.CPL)
##              Df  Sum Sq  Mean Sq F value   Pr(>F)    
## BASIN         5 0.02013 0.004025   5.383 0.000159 ***
## Residuals   128 0.09572 0.000748                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
A1.PreDL <- aov(abs.res.PreDL ~ BASIN, data=form3)
summary(A1.PreDL)
##              Df Sum Sq  Mean Sq F value Pr(>F)
## BASIN         5 0.0222 0.004437   1.299  0.268
## Residuals   128 0.4372 0.003415
ggplot(form3, aes(x=factor(BASIN), y=abs.res.PreDL)) +
  geom_boxplot()

A1.DbL <- aov(abs.res.DbL ~ BASIN, data=form3)
summary(A1.DbL)
##              Df  Sum Sq  Mean Sq F value   Pr(>F)    
## BASIN         5 0.03427 0.006855   8.464 6.05e-07 ***
## Residuals   128 0.10367 0.000810                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ggplot(form3, aes(x=factor(BASIN), y=abs.res.DbL)) +
  geom_boxplot()

A1.HL <- aov(abs.res.HL ~ BASIN, data=form3)
summary(A1.HL)
##              Df  Sum Sq   Mean Sq F value Pr(>F)
## BASIN         5 0.00766 0.0015321   1.809  0.115
## Residuals   128 0.10838 0.0008467
ggplot(form3, aes(x=factor(BASIN), y=abs.res.HL)) +
  geom_boxplot()

A1.HD <- aov(abs.res.HD ~ BASIN, data=form3)
summary(A1.HD)
##              Df  Sum Sq  Mean Sq F value Pr(>F)
## BASIN         5 0.00684 0.001367   1.228    0.3
## Residuals   128 0.14257 0.001114
ggplot(form3, aes(x=factor(BASIN), y=abs.res.HD)) +
  geom_boxplot()

A1.HW <- aov(abs.res.HW ~ BASIN, data=form3)
summary(A1.HW)
##              Df  Sum Sq   Mean Sq F value   Pr(>F)    
## BASIN         5 0.01460 0.0029194   5.746 8.11e-05 ***
## Residuals   128 0.06504 0.0005081                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ggplot(form3, aes(x=factor(BASIN), y=abs.res.HW)) +
  geom_boxplot()

A1.SnL <- aov(abs.res.SnL ~ BASIN, data=form3)
summary(A1.SnL)
##              Df  Sum Sq  Mean Sq F value   Pr(>F)    
## BASIN         5 0.01931 0.003861   9.138 1.87e-07 ***
## Residuals   128 0.05409 0.000423                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
A1.OL <- aov(abs.res.OL ~ BASIN, data=form3)
summary(A1.OL)
##              Df    Sum Sq   Mean Sq F value Pr(>F)
## BASIN         5 1.240e-32 2.479e-33   0.511  0.768
## Residuals   128 6.214e-31 4.854e-33
ggplot(form3, aes(x=factor(BASIN), y=abs.res.OL)) +
  geom_boxplot()

Watersheds
A2.D <- aov(abs.res.D ~ WATERSHED, data=form3)
summary(A2.D)
##              Df Sum Sq Mean Sq F value Pr(>F)  
## WATERSHED    11  3.847  0.3497   1.963 0.0378 *
## Residuals   122 21.736  0.1782                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ggplot(form3, aes(x=factor(WATERSHED), y=abs.res.D)) +
  geom_boxplot()

A2.P1 <- aov(abs.res.P1 ~ WATERSHED, data=form3)
summary(A2.P1)
##              Df Sum Sq Mean Sq F value Pr(>F)
## WATERSHED    11  1.302  0.1184   1.005  0.446
## Residuals   122 14.368  0.1178
A2.P1.R <- aov(abs.res.P1.R ~ WATERSHED, data=form3)
summary(A2.P1.R)
##              Df Sum Sq Mean Sq F value Pr(>F)
## WATERSHED    11   0.34 0.03093   0.226  0.995
## Residuals   122  16.71 0.13698
ggplot(form3, aes(x=factor(WATERSHED), y=abs.res.P1.R)) +
  geom_boxplot()

A2.LLSC <- aov(abs.res.LLSC ~ WATERSHED, data=form3)
summary(A2.LLSC)
##              Df Sum Sq Mean Sq F value Pr(>F)  
## WATERSHED    11   8.16  0.7422   2.267  0.015 *
## Residuals   122  39.94  0.3274                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ggplot(form3, aes(x=factor(WATERSHED), y=abs.res.LLSC)) +
  geom_boxplot()

A2.SBLL <- aov(abs.res.SBLL ~ WATERSHED, data=form3)
summary(A2.SBLL)
##              Df Sum Sq Mean Sq F value   Pr(>F)    
## WATERSHED    11  5.086  0.4623   4.741 4.93e-06 ***
## Residuals   122 11.897  0.0975                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ggplot(form3, aes(x=factor(WATERSHED), y=abs.res.SBLL)) +
  geom_boxplot()

A2.BD <- aov(abs.res.BD ~ WATERSHED, data=form3)
summary(A2.BD)
##              Df  Sum Sq   Mean Sq F value Pr(>F)
## WATERSHED    11 0.00851 0.0007737   0.944  0.501
## Residuals   122 0.09999 0.0008196
ggplot(form3, aes(x=factor(WATERSHED), y=abs.res.BD)) +
  geom_boxplot()

A2.CPD <- aov(abs.res.CPD ~ WATERSHED, data=form3)
summary(A2.CPD)
##              Df  Sum Sq   Mean Sq F value  Pr(>F)   
## WATERSHED    11 0.01833 0.0016662    2.77 0.00304 **
## Residuals   122 0.07338 0.0006015                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ggplot(form3, aes(x=factor(WATERSHED), y=abs.res.CPD)) +
  geom_boxplot()

A2.CPL <- aov(abs.res.CPL ~ WATERSHED, data=form3)
summary(A2.CPL)
##              Df  Sum Sq   Mean Sq F value   Pr(>F)    
## WATERSHED    11 0.03237 0.0029429   4.301 2.05e-05 ***
## Residuals   122 0.08347 0.0006842                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
A2.PreDL <- aov(abs.res.PreDL ~ WATERSHED, data=form3)
summary(A2.PreDL)
##              Df Sum Sq  Mean Sq F value Pr(>F)
## WATERSHED    11 0.0291 0.002646    0.75  0.688
## Residuals   122 0.4303 0.003527
A2.DbL <- aov(abs.res.DbL ~ WATERSHED, data=form3)
summary(A2.DbL)
##              Df  Sum Sq  Mean Sq F value   Pr(>F)    
## WATERSHED    11 0.03993 0.003630   4.518 1.01e-05 ***
## Residuals   122 0.09801 0.000803                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ggplot(form3, aes(x=factor(WATERSHED), y=abs.res.DbL)) +
  geom_boxplot()

A2.HL <- aov(abs.res.HL ~ WATERSHED, data=form3)
summary(A2.HL)
##              Df  Sum Sq   Mean Sq F value Pr(>F)
## WATERSHED    11 0.01169 0.0010627   1.243  0.267
## Residuals   122 0.10435 0.0008553
ggplot(form3, aes(x=factor(WATERSHED), y=abs.res.HL)) +
  geom_boxplot()

A2.HD <- aov(abs.res.HD ~ WATERSHED, data=form3)
summary(A2.HD)
##              Df  Sum Sq   Mean Sq F value Pr(>F)
## WATERSHED    11 0.01036 0.0009418   0.826  0.614
## Residuals   122 0.13905 0.0011397
ggplot(form3, aes(x=factor(WATERSHED), y=abs.res.HD)) +
  geom_boxplot()

A2.HW <- aov(abs.res.HW ~ WATERSHED, data=form3)
summary(A2.HW)
##              Df  Sum Sq   Mean Sq F value  Pr(>F)   
## WATERSHED    11 0.01704 0.0015488   3.018 0.00136 **
## Residuals   122 0.06260 0.0005131                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ggplot(form3, aes(x=factor(WATERSHED), y=abs.res.HW)) +
  geom_boxplot()

A2.SnL <- aov(abs.res.SnL ~ WATERSHED, data=form3)
summary(A2.SnL)
##              Df  Sum Sq   Mean Sq F value   Pr(>F)    
## WATERSHED    11 0.02208 0.0020073   4.772 4.46e-06 ***
## Residuals   122 0.05131 0.0004206                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
A2.OL <- aov(abs.res.OL ~ WATERSHED, data=form3)
summary(A2.OL)
##              Df    Sum Sq   Mean Sq F value Pr(>F)
## WATERSHED    11 2.610e-32 2.372e-33   0.476  0.915
## Residuals   122 6.077e-31 4.981e-33
ggplot(form3, aes(x=factor(WATERSHED), y=abs.res.OL)) +
  geom_boxplot()

Sailfin

Zones
library(ggplot2)

A.D <- aov(abs.res.D ~ QUARTILE, data=lat3)
summary(A.D)
##              Df Sum Sq Mean Sq F value  Pr(>F)   
## QUARTILE      3  2.306  0.7685   4.292 0.00635 **
## Residuals   130 23.277  0.1791                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ggplot(lat3, aes(x=factor(QUARTILE), y=abs.res.D)) +
  geom_boxplot()

A.P1 <- aov(abs.res.P1 ~ QUARTILE, data=lat3)
summary(A.P1)
##              Df Sum Sq Mean Sq F value Pr(>F)
## QUARTILE      3  0.159 0.05296   0.444  0.722
## Residuals   130 15.511 0.11932
A.P1.R <- aov(abs.res.P1.R ~ QUARTILE, data=lat3)
summary(A.P1.R)
##              Df Sum Sq Mean Sq F value Pr(>F)
## QUARTILE      3  0.065 0.02172   0.166  0.919
## Residuals   130 16.987 0.13067
ggplot(lat3, aes(x=factor(QUARTILE), y=abs.res.P1.R)) +
  geom_boxplot()

A.LLSC <- aov(abs.res.LLSC ~ QUARTILE, data=lat3)
summary(A.LLSC)
##              Df Sum Sq Mean Sq F value Pr(>F)
## QUARTILE      3   1.31  0.4354   1.209  0.309
## Residuals   130  46.80  0.3600
A.SBLL <- aov(abs.res.SBLL ~ QUARTILE, data=lat3)
summary(A.SBLL)
##              Df Sum Sq Mean Sq F value Pr(>F)
## QUARTILE      3  0.436  0.1454   1.142  0.335
## Residuals   130 16.546  0.1273
ggplot(lat3, aes(x=factor(QUARTILE), y=abs.res.SBLL)) +
  geom_boxplot()

A.BD <- aov(abs.res.BD ~ QUARTILE, data=lat3)
summary(A.BD)
##              Df  Sum Sq   Mean Sq F value Pr(>F)
## QUARTILE      3 0.00259 0.0008622   1.058  0.369
## Residuals   130 0.10591 0.0008147
ggplot(lat3, aes(x=factor(QUARTILE), y=abs.res.BD)) +
  geom_boxplot()

A.CPD <- aov(abs.res.CPD ~ QUARTILE, data=lat3)
summary(A.CPD)
##              Df  Sum Sq  Mean Sq F value  Pr(>F)    
## QUARTILE      3 0.01249 0.004163   6.832 0.00026 ***
## Residuals   130 0.07922 0.000609                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ggplot(lat3, aes(x=factor(QUARTILE), y=abs.res.CPD)) +
  geom_boxplot()

A.CPL <- aov(abs.res.CPL ~ QUARTILE, data=lat3)
summary(A.CPL)
##              Df  Sum Sq  Mean Sq F value   Pr(>F)    
## QUARTILE      3 0.02074 0.006912   9.448 1.08e-05 ***
## Residuals   130 0.09511 0.000732                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ggplot(lat3, aes(x=factor(QUARTILE), y=abs.res.CPL)) +
  geom_boxplot()

A.PreDL <- aov(abs.res.PreDL ~ QUARTILE, data=lat3)
summary(A.PreDL)
##              Df Sum Sq  Mean Sq F value Pr(>F)
## QUARTILE      3 0.0082 0.002742    0.79  0.502
## Residuals   130 0.4511 0.003470
ggplot(lat3, aes(x=factor(QUARTILE), y=abs.res.PreDL)) +
  geom_boxplot()

A.DbL <- aov(abs.res.DbL ~ QUARTILE, data=lat3)
summary(A.DbL)
##              Df Sum Sq  Mean Sq F value   Pr(>F)    
## QUARTILE      3 0.0195 0.006501   7.136 0.000179 ***
## Residuals   130 0.1184 0.000911                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ggplot(lat3, aes(x=factor(QUARTILE), y=abs.res.DbL)) +
  geom_boxplot()

A.HL <- aov(abs.res.HL ~ QUARTILE, data=lat3)
summary(A.HL)
##              Df  Sum Sq   Mean Sq F value Pr(>F)
## QUARTILE      3 0.00344 0.0011466   1.324  0.269
## Residuals   130 0.11260 0.0008661
ggplot(lat3, aes(x=factor(QUARTILE), y=abs.res.HL)) +
  geom_boxplot()

A.HD <- aov(abs.res.HD ~ QUARTILE, data=lat3)
summary(A.HD)
##              Df  Sum Sq  Mean Sq F value Pr(>F)
## QUARTILE      3 0.00672 0.002239    2.04  0.111
## Residuals   130 0.14269 0.001098
A.HW <- aov(abs.res.HW ~ QUARTILE, data=lat3)
summary(A.HW)
##              Df  Sum Sq  Mean Sq F value   Pr(>F)    
## QUARTILE      3 0.01136 0.003788   7.212 0.000163 ***
## Residuals   130 0.06827 0.000525                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ggplot(lat3, aes(x=factor(QUARTILE), y=abs.res.HW)) +
  geom_boxplot()

A.SnL <- aov(abs.res.SnL ~ QUARTILE, data=lat3)
summary(A.SnL)
##              Df  Sum Sq   Mean Sq F value  Pr(>F)   
## QUARTILE      3 0.00766 0.0025545   5.052 0.00242 **
## Residuals   130 0.06573 0.0005056                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
A.OL <- aov(abs.res.OL ~ QUARTILE, data=lat3)
summary(A.OL)
##              Df   Sum Sq   Mean Sq F value Pr(>F)
## QUARTILE      3 2.17e-32 7.242e-33   1.538  0.208
## Residuals   130 6.12e-31 4.708e-33
ggplot(lat3, aes(x=factor(QUARTILE), y=abs.res.OL)) +
  geom_boxplot()

Basins
A1.D <- aov(abs.res.D ~ BASIN, data=lat3)
summary(A1.D)
##              Df Sum Sq Mean Sq F value Pr(>F)  
## BASIN         5  2.789  0.5579   3.133 0.0106 *
## Residuals   128 22.793  0.1781                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ggplot(lat3, aes(x=factor(BASIN), y=abs.res.D)) +
  geom_boxplot()

A1.P1 <- aov(abs.res.P1 ~ BASIN, data=lat3)
summary(A1.P1)
##              Df Sum Sq Mean Sq F value Pr(>F)
## BASIN         5  0.945  0.1891   1.643  0.153
## Residuals   128 14.725  0.1150
A1.P1.R <- aov(abs.res.P1.R ~ BASIN, data=lat3)
summary(A1.P1.R)
##              Df Sum Sq Mean Sq F value Pr(>F)
## BASIN         5  0.159 0.03187   0.241  0.943
## Residuals   128 16.893 0.13197
ggplot(lat3, aes(x=factor(BASIN), y=abs.res.P1.R)) +
  geom_boxplot()

A1.LLSC <- aov(abs.res.LLSC ~ BASIN, data=lat3)
summary(A1.LLSC)
##              Df Sum Sq Mean Sq F value Pr(>F)  
## BASIN         5   4.54  0.9080   2.668  0.025 *
## Residuals   128  43.57  0.3404                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ggplot(lat3, aes(x=factor(BASIN), y=abs.res.LLSC)) +
  geom_boxplot()

A1.SBLL <- aov(abs.res.SBLL ~ BASIN, data=lat3)
summary(A1.SBLL)
##              Df Sum Sq Mean Sq F value Pr(>F)
## BASIN         5  0.755  0.1510   1.191  0.317
## Residuals   128 16.228  0.1268
ggplot(lat3, aes(x=factor(BASIN), y=abs.res.SBLL)) +
  geom_boxplot()

A1.BD <- aov(abs.res.BD ~ BASIN, data=lat3)
summary(A1.BD)
##              Df  Sum Sq   Mean Sq F value Pr(>F)
## BASIN         5 0.00316 0.0006313   0.767  0.575
## Residuals   128 0.10534 0.0008230
ggplot(lat3, aes(x=factor(BASIN), y=abs.res.BD)) +
  geom_boxplot()

A1.CPD <- aov(abs.res.CPD ~ BASIN, data=lat3)
summary(A1.CPD)
##              Df  Sum Sq  Mean Sq F value   Pr(>F)    
## BASIN         5 0.01713 0.003425   5.879 6.34e-05 ***
## Residuals   128 0.07458 0.000583                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ggplot(lat3, aes(x=factor(BASIN), y=abs.res.CPD)) +
  geom_boxplot()

A1.CPL <- aov(abs.res.CPL ~ BASIN, data=lat3)
summary(A1.CPL)
##              Df  Sum Sq  Mean Sq F value   Pr(>F)    
## BASIN         5 0.02013 0.004025   5.383 0.000159 ***
## Residuals   128 0.09572 0.000748                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
A1.PreDL <- aov(abs.res.PreDL ~ BASIN, data=lat3)
summary(A1.PreDL)
##              Df Sum Sq  Mean Sq F value Pr(>F)
## BASIN         5 0.0222 0.004437   1.299  0.268
## Residuals   128 0.4372 0.003415
ggplot(lat3, aes(x=factor(BASIN), y=abs.res.PreDL)) +
  geom_boxplot()

A1.DbL <- aov(abs.res.DbL ~ BASIN, data=lat3)
summary(A1.DbL)
##              Df  Sum Sq  Mean Sq F value   Pr(>F)    
## BASIN         5 0.03427 0.006855   8.464 6.05e-07 ***
## Residuals   128 0.10367 0.000810                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ggplot(lat3, aes(x=factor(BASIN), y=abs.res.DbL)) +
  geom_boxplot()

A1.HL <- aov(abs.res.HL ~ BASIN, data=lat3)
summary(A1.HL)
##              Df  Sum Sq   Mean Sq F value Pr(>F)
## BASIN         5 0.00766 0.0015321   1.809  0.115
## Residuals   128 0.10838 0.0008467
ggplot(lat3, aes(x=factor(BASIN), y=abs.res.HL)) +
  geom_boxplot()

A1.HD <- aov(abs.res.HD ~ BASIN, data=lat3)
summary(A1.HD)
##              Df  Sum Sq  Mean Sq F value Pr(>F)
## BASIN         5 0.00684 0.001367   1.228    0.3
## Residuals   128 0.14257 0.001114
ggplot(lat3, aes(x=factor(BASIN), y=abs.res.HD)) +
  geom_boxplot()

A1.HW <- aov(abs.res.HW ~ BASIN, data=lat3)
summary(A1.HW)
##              Df  Sum Sq   Mean Sq F value   Pr(>F)    
## BASIN         5 0.01460 0.0029194   5.746 8.11e-05 ***
## Residuals   128 0.06504 0.0005081                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ggplot(lat3, aes(x=factor(BASIN), y=abs.res.HW)) +
  geom_boxplot()

A1.SnL <- aov(abs.res.SnL ~ BASIN, data=lat3)
summary(A1.SnL)
##              Df  Sum Sq  Mean Sq F value   Pr(>F)    
## BASIN         5 0.01931 0.003861   9.138 1.87e-07 ***
## Residuals   128 0.05409 0.000423                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
A1.OL <- aov(abs.res.OL ~ BASIN, data=lat3)
summary(A1.OL)
##              Df    Sum Sq   Mean Sq F value Pr(>F)
## BASIN         5 1.240e-32 2.479e-33   0.511  0.768
## Residuals   128 6.214e-31 4.854e-33
ggplot(lat3, aes(x=factor(BASIN), y=abs.res.OL)) +
  geom_boxplot()

Watersheds
A2.D <- aov(abs.res.D ~ WATERSHED, data=lat3)
summary(A2.D)
##              Df Sum Sq Mean Sq F value Pr(>F)  
## WATERSHED    11  3.847  0.3497   1.963 0.0378 *
## Residuals   122 21.736  0.1782                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ggplot(lat3, aes(x=factor(WATERSHED), y=abs.res.D)) +
  geom_boxplot()

A2.P1 <- aov(abs.res.P1 ~ WATERSHED, data=lat3)
summary(A2.P1)
##              Df Sum Sq Mean Sq F value Pr(>F)
## WATERSHED    11  1.302  0.1184   1.005  0.446
## Residuals   122 14.368  0.1178
A2.P1.R <- aov(abs.res.P1.R ~ WATERSHED, data=lat3)
summary(A2.P1.R)
##              Df Sum Sq Mean Sq F value Pr(>F)
## WATERSHED    11   0.34 0.03093   0.226  0.995
## Residuals   122  16.71 0.13698
ggplot(lat3, aes(x=factor(WATERSHED), y=abs.res.P1.R)) +
  geom_boxplot()

A2.LLSC <- aov(abs.res.LLSC ~ WATERSHED, data=lat3)
summary(A2.LLSC)
##              Df Sum Sq Mean Sq F value Pr(>F)  
## WATERSHED    11   8.16  0.7422   2.267  0.015 *
## Residuals   122  39.94  0.3274                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ggplot(lat3, aes(x=factor(WATERSHED), y=abs.res.LLSC)) +
  geom_boxplot()

A2.SBLL <- aov(abs.res.SBLL ~ WATERSHED, data=lat3)
summary(A2.SBLL)
##              Df Sum Sq Mean Sq F value   Pr(>F)    
## WATERSHED    11  5.086  0.4623   4.741 4.93e-06 ***
## Residuals   122 11.897  0.0975                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ggplot(lat3, aes(x=factor(WATERSHED), y=abs.res.SBLL)) +
  geom_boxplot()

A2.BD <- aov(abs.res.BD ~ WATERSHED, data=lat3)
summary(A2.BD)
##              Df  Sum Sq   Mean Sq F value Pr(>F)
## WATERSHED    11 0.00851 0.0007737   0.944  0.501
## Residuals   122 0.09999 0.0008196
ggplot(lat3, aes(x=factor(WATERSHED), y=abs.res.BD)) +
  geom_boxplot()

A2.CPD <- aov(abs.res.CPD ~ WATERSHED, data=lat3)
summary(A2.CPD)
##              Df  Sum Sq   Mean Sq F value  Pr(>F)   
## WATERSHED    11 0.01833 0.0016662    2.77 0.00304 **
## Residuals   122 0.07338 0.0006015                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ggplot(lat3, aes(x=factor(WATERSHED), y=abs.res.CPD)) +
  geom_boxplot()

A2.CPL <- aov(abs.res.CPL ~ WATERSHED, data=lat3)
summary(A2.CPL)
##              Df  Sum Sq   Mean Sq F value   Pr(>F)    
## WATERSHED    11 0.03237 0.0029429   4.301 2.05e-05 ***
## Residuals   122 0.08347 0.0006842                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
A2.PreDL <- aov(abs.res.PreDL ~ WATERSHED, data=lat3)
summary(A2.PreDL)
##              Df Sum Sq  Mean Sq F value Pr(>F)
## WATERSHED    11 0.0291 0.002646    0.75  0.688
## Residuals   122 0.4303 0.003527
A2.DbL <- aov(abs.res.DbL ~ WATERSHED, data=lat3)
summary(A2.DbL)
##              Df  Sum Sq  Mean Sq F value   Pr(>F)    
## WATERSHED    11 0.03993 0.003630   4.518 1.01e-05 ***
## Residuals   122 0.09801 0.000803                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ggplot(lat3, aes(x=factor(WATERSHED), y=abs.res.DbL)) +
  geom_boxplot()

A2.HL <- aov(abs.res.HL ~ WATERSHED, data=lat3)
summary(A2.HL)
##              Df  Sum Sq   Mean Sq F value Pr(>F)
## WATERSHED    11 0.01169 0.0010627   1.243  0.267
## Residuals   122 0.10435 0.0008553
ggplot(lat3, aes(x=factor(WATERSHED), y=abs.res.HL)) +
  geom_boxplot()

A2.HD <- aov(abs.res.HD ~ WATERSHED, data=lat3)
summary(A2.HD)
##              Df  Sum Sq   Mean Sq F value Pr(>F)
## WATERSHED    11 0.01036 0.0009418   0.826  0.614
## Residuals   122 0.13905 0.0011397
ggplot(lat3, aes(x=factor(WATERSHED), y=abs.res.HD)) +
  geom_boxplot()

A2.HW <- aov(abs.res.HW ~ WATERSHED, data=lat3)
summary(A2.HW)
##              Df  Sum Sq   Mean Sq F value  Pr(>F)   
## WATERSHED    11 0.01704 0.0015488   3.018 0.00136 **
## Residuals   122 0.06260 0.0005131                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ggplot(lat3, aes(x=factor(WATERSHED), y=abs.res.HW)) +
  geom_boxplot()

A2.SnL <- aov(abs.res.SnL ~ WATERSHED, data=lat3)
summary(A2.SnL)
##              Df  Sum Sq   Mean Sq F value   Pr(>F)    
## WATERSHED    11 0.02208 0.0020073   4.772 4.46e-06 ***
## Residuals   122 0.05131 0.0004206                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
A2.OL <- aov(abs.res.OL ~ WATERSHED, data=lat3)
summary(A2.OL)
##              Df    Sum Sq   Mean Sq F value Pr(>F)
## WATERSHED    11 2.610e-32 2.372e-33   0.476  0.915
## Residuals   122 6.077e-31 4.981e-33
ggplot(lat3, aes(x=factor(WATERSHED), y=abs.res.OL)) +
  geom_boxplot()

Both

Zones
library(ggplot2)

A.D <- aov(abs.res.D ~ SPP*QUARTILE, data=raw3)
summary(A.D)
##               Df Sum Sq Mean Sq F value  Pr(>F)    
## SPP            1   0.06  0.0635   0.408 0.52344    
## QUARTILE       3   0.45  0.1510   0.970 0.40715    
## SPP:QUARTILE   3   3.22  1.0719   6.888 0.00017 ***
## Residuals    292  45.44  0.1556                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ggplot(raw3, aes(x=factor(QUARTILE), y=abs.res.D, fill=SPP)) +
  geom_boxplot()

A.P1 <- aov(abs.res.P1 ~ SPP*QUARTILE, data=raw3)
summary(A.P1)
##               Df Sum Sq Mean Sq F value Pr(>F)
## SPP            1   0.29  0.2878   1.554  0.214
## QUARTILE       3   1.10  0.3663   1.977  0.117
## SPP:QUARTILE   3   0.73  0.2427   1.310  0.271
## Residuals    292  54.10  0.1853
A.P1.R <- aov(abs.res.P1.R ~ SPP*QUARTILE, data=raw3)
summary(A.P1.R)
##               Df Sum Sq Mean Sq F value   Pr(>F)    
## SPP            1   1.80  1.7983  11.423 0.000824 ***
## QUARTILE       3   0.13  0.0434   0.276 0.842918    
## SPP:QUARTILE   3   0.06  0.0189   0.120 0.948162    
## Residuals    292  45.97  0.1574                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ggplot(raw3, aes(x=factor(QUARTILE), y=abs.res.P1.R, fill=SPP)) +
  geom_boxplot()

A.LLSC <- aov(abs.res.LLSC ~ SPP*QUARTILE, data=raw3)
summary(A.LLSC)
##               Df Sum Sq Mean Sq F value Pr(>F)  
## SPP            1   1.40  1.4013   2.940 0.0875 .
## QUARTILE       3   1.35  0.4506   0.945 0.4190  
## SPP:QUARTILE   3   3.42  1.1399   2.391 0.0688 .
## Residuals    292 139.19  0.4767                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
A.SBLL <- aov(abs.res.SBLL ~ SPP*QUARTILE, data=raw3)
summary(A.SBLL)
##               Df Sum Sq Mean Sq F value Pr(>F)  
## SPP            1   0.65  0.6523   4.859 0.0283 *
## QUARTILE       3   0.78  0.2591   1.930 0.1248  
## SPP:QUARTILE   3   0.12  0.0390   0.291 0.8320  
## Residuals    292  39.20  0.1342                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ggplot(raw3, aes(x=factor(QUARTILE), y=abs.res.SBLL, fill=SPP)) +
  geom_boxplot()

A.BD <- aov(abs.res.BD ~ SPP*QUARTILE, data=raw3)
summary(A.BD)
##               Df Sum Sq  Mean Sq F value Pr(>F)  
## SPP            1 0.0023 0.002282   0.608 0.4363  
## QUARTILE       3 0.0161 0.005364   1.429 0.2345  
## SPP:QUARTILE   3 0.0289 0.009643   2.568 0.0546 .
## Residuals    292 1.0963 0.003755                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ggplot(raw3, aes(x=factor(QUARTILE), y=abs.res.BD, fill=SPP)) +
  geom_boxplot()

A.CPD <- aov(abs.res.CPD ~ SPP*QUARTILE, data=raw3)
summary(A.CPD)
##               Df  Sum Sq   Mean Sq F value  Pr(>F)   
## SPP            1 0.00312 0.0031247   4.900 0.02764 * 
## QUARTILE       3 0.00786 0.0026201   4.108 0.00707 **
## SPP:QUARTILE   3 0.00713 0.0023775   3.728 0.01176 * 
## Residuals    292 0.18622 0.0006377                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ggplot(raw3, aes(x=factor(QUARTILE), y=abs.res.CPD, fill=SPP)) +
  geom_boxplot()

A.CPL <- aov(abs.res.CPL ~ SPP*QUARTILE, data=raw3)
summary(A.CPL)
##               Df  Sum Sq  Mean Sq F value   Pr(>F)    
## SPP            1 0.00000 0.000000   0.000 0.986836    
## QUARTILE       3 0.01640 0.005468   6.588 0.000254 ***
## SPP:QUARTILE   3 0.02033 0.006777   8.166 3.08e-05 ***
## Residuals    292 0.24236 0.000830                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ggplot(raw3, aes(x=factor(QUARTILE), y=abs.res.CPL, fill=SPP)) +
  geom_boxplot()

A.PreDL <- aov(abs.res.PreDL ~ SPP*QUARTILE, data=raw3)
summary(A.PreDL)
##               Df Sum Sq  Mean Sq F value Pr(>F)
## SPP            1 0.0012 0.001166   0.501  0.480
## QUARTILE       3 0.0053 0.001774   0.763  0.516
## SPP:QUARTILE   3 0.0055 0.001820   0.783  0.504
## Residuals    292 0.6792 0.002326
ggplot(raw3, aes(x=factor(QUARTILE), y=abs.res.PreDL, fill=SPP)) +
  geom_boxplot()

A.DbL <- aov(abs.res.DbL ~ SPP*QUARTILE, data=raw3)
summary(A.DbL)
##               Df  Sum Sq  Mean Sq F value  Pr(>F)   
## SPP            1 0.00556 0.005557   5.888 0.01585 * 
## QUARTILE       3 0.01021 0.003404   3.607 0.01382 * 
## SPP:QUARTILE   3 0.01319 0.004397   4.659 0.00338 **
## Residuals    292 0.27557 0.000944                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ggplot(raw3, aes(x=factor(QUARTILE), y=abs.res.DbL, fill=SPP)) +
  geom_boxplot()

A.HL <- aov(abs.res.HL ~ SPP*QUARTILE, data=raw3)
summary(A.HL)
##               Df  Sum Sq  Mean Sq F value   Pr(>F)    
## SPP            1 0.00081 0.000812   0.861   0.3541    
## QUARTILE       3 0.02099 0.006995   7.421 8.33e-05 ***
## SPP:QUARTILE   3 0.01029 0.003432   3.640   0.0132 *  
## Residuals    292 0.27527 0.000943                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ggplot(raw3, aes(x=factor(QUARTILE), y=abs.res.HL, fill=SPP)) +
  geom_boxplot()

A.HD <- aov(abs.res.HD ~ SPP*QUARTILE, data=raw3)
summary(A.HD)
##               Df  Sum Sq   Mean Sq F value Pr(>F)
## SPP            1 0.00000 0.0000000   0.000  0.999
## QUARTILE       3 0.00443 0.0014757   1.522  0.209
## SPP:QUARTILE   3 0.00582 0.0019390   1.999  0.114
## Residuals    292 0.28319 0.0009698
A.HW <- aov(abs.res.HW ~ SPP*QUARTILE, data=raw3)
summary(A.HW)
##               Df  Sum Sq   Mean Sq F value  Pr(>F)   
## SPP            1 0.00213 0.0021307   2.861 0.09184 . 
## QUARTILE       3 0.00917 0.0030561   4.103 0.00712 **
## SPP:QUARTILE   3 0.00710 0.0023664   3.177 0.02448 * 
## Residuals    292 0.21750 0.0007449                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ggplot(raw3, aes(x=factor(QUARTILE), y=abs.res.HW, fill=SPP)) +
  geom_boxplot()

A.SnL <- aov(abs.res.SnL ~ SPP*QUARTILE, data=raw3)
summary(A.SnL)
##               Df  Sum Sq  Mean Sq F value   Pr(>F)    
## SPP            1 0.00009 0.000095   0.163   0.6867    
## QUARTILE       3 0.01743 0.005809   9.995 2.73e-06 ***
## SPP:QUARTILE   3 0.00494 0.001647   2.833   0.0386 *  
## Residuals    292 0.16972 0.000581                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
A.OL <- aov(abs.res.OL ~ SPP*QUARTILE, data=raw3)
summary(A.OL)
##               Df    Sum Sq   Mean Sq F value Pr(>F)  
## SPP            1 4.900e-33 4.915e-33   2.108 0.1476  
## QUARTILE       3 1.830e-32 6.103e-33   2.617 0.0512 .
## SPP:QUARTILE   3 5.800e-33 1.948e-33   0.835 0.4754  
## Residuals    292 6.809e-31 2.332e-33                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ggplot(raw3, aes(x=factor(QUARTILE), y=abs.res.OL, fill=SPP)) +
  geom_boxplot()

Basins
A1.D <- aov(abs.res.D ~ SPP*BASIN, data=raw3)
summary(A1.D)
##              Df Sum Sq Mean Sq F value  Pr(>F)   
## SPP           1   0.06  0.0635   0.402 0.52666   
## BASIN         6   3.12  0.5205   3.293 0.00377 **
## SPP:BASIN     3   0.31  0.1028   0.650 0.58344   
## Residuals   289  45.68  0.1581                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ggplot(raw3, aes(x=factor(BASIN), y=abs.res.D, fill=SPP)) +
  geom_boxplot()

A1.P1 <- aov(abs.res.P1 ~ SPP*BASIN, data=raw3)
summary(A1.P1)
##              Df Sum Sq Mean Sq F value Pr(>F)
## SPP           1   0.29  0.2878   1.528  0.217
## BASIN         6   0.97  0.1613   0.856  0.528
## SPP:BASIN     3   0.53  0.1751   0.930  0.427
## Residuals   289  54.43  0.1883
A1.P1.R <- aov(abs.res.P1.R ~ SPP*BASIN, data=raw3)
summary(A1.P1.R)
##              Df Sum Sq Mean Sq F value   Pr(>F)    
## SPP           1   1.80  1.7983  11.421 0.000826 ***
## BASIN         6   0.51  0.0845   0.537 0.780341    
## SPP:BASIN     3   0.15  0.0485   0.308 0.819718    
## Residuals   289  45.50  0.1575                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ggplot(raw3, aes(x=factor(BASIN), y=abs.res.P1.R, fill=SPP)) +
  geom_boxplot()

A1.LLSC <- aov(abs.res.LLSC ~ SPP*BASIN, data=raw3)
summary(A1.LLSC)
##              Df Sum Sq Mean Sq F value Pr(>F)  
## SPP           1   1.40  1.4013   2.999 0.0844 .
## BASIN         6   6.83  1.1376   2.435 0.0260 *
## SPP:BASIN     3   2.11  0.7026   1.504 0.2137  
## Residuals   289 135.02  0.4672                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ggplot(raw3, aes(x=factor(BASIN), y=abs.res.LLSC, fill=SPP)) +
  geom_boxplot()

A1.SBLL <- aov(abs.res.SBLL ~ SPP*BASIN, data=raw3)
summary(A1.SBLL)
##              Df Sum Sq Mean Sq F value Pr(>F)  
## SPP           1   0.65  0.6523   4.895 0.0277 *
## BASIN         6   1.07  0.1785   1.340 0.2393  
## SPP:BASIN     3   0.51  0.1716   1.288 0.2788  
## Residuals   289  38.51  0.1332                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ggplot(raw3, aes(x=factor(BASIN), y=abs.res.SBLL, fill=SPP)) +
  geom_boxplot()

A1.BD <- aov(abs.res.BD ~ SPP*BASIN, data=raw3)
summary(A1.BD)
##              Df Sum Sq  Mean Sq F value Pr(>F)
## SPP           1 0.0023 0.002282   0.590  0.443
## BASIN         6 0.0161 0.002675   0.691  0.657
## SPP:BASIN     3 0.0068 0.002270   0.586  0.624
## Residuals   289 1.1185 0.003870
ggplot(raw3, aes(x=factor(BASIN), y=abs.res.BD, fill=SPP)) +
  geom_boxplot()

A1.CPD <- aov(abs.res.CPD ~ SPP*BASIN, data=raw3)
summary(A1.CPD)
##              Df  Sum Sq   Mean Sq F value   Pr(>F)    
## SPP           1 0.00312 0.0031247   4.984 0.026343 *  
## BASIN         6 0.01576 0.0026273   4.191 0.000464 ***
## SPP:BASIN     3 0.00428 0.0014261   2.275 0.080070 .  
## Residuals   289 0.18117 0.0006269                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ggplot(raw3, aes(x=factor(BASIN), y=abs.res.CPD, fill=SPP)) +
  geom_boxplot()

A1.CPL <- aov(abs.res.CPL ~ SPP*BASIN, data=raw3)
summary(A1.CPL)
##              Df  Sum Sq  Mean Sq F value   Pr(>F)    
## SPP           1 0.00000 0.000000   0.000    0.987    
## BASIN         6 0.03228 0.005380   6.320 2.92e-06 ***
## SPP:BASIN     3 0.00078 0.000259   0.304    0.822    
## Residuals   289 0.24603 0.000851                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
A1.PreDL <- aov(abs.res.PreDL ~ SPP*BASIN, data=raw3)
summary(A1.PreDL)
##              Df Sum Sq  Mean Sq F value Pr(>F)
## SPP           1 0.0012 0.001166   0.507  0.477
## BASIN         6 0.0132 0.002195   0.954  0.456
## SPP:BASIN     3 0.0122 0.004078   1.773  0.152
## Residuals   289 0.6646 0.002300
ggplot(raw3, aes(x=factor(BASIN), y=abs.res.PreDL, fill=SPP)) +
  geom_boxplot()

A1.DbL <- aov(abs.res.DbL ~ SPP*BASIN, data=raw3)
summary(A1.DbL)
##              Df  Sum Sq  Mean Sq F value   Pr(>F)    
## SPP           1 0.00556 0.005557   6.191 0.013405 *  
## BASIN         6 0.02192 0.003653   4.070 0.000618 ***
## SPP:BASIN     3 0.01766 0.005886   6.557 0.000265 ***
## Residuals   289 0.25940 0.000898                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ggplot(raw3, aes(x=factor(BASIN), y=abs.res.DbL, fill=SPP)) +
  geom_boxplot()

A1.HL <- aov(abs.res.HL ~ SPP*BASIN, data=raw3)
summary(A1.HL)
##              Df  Sum Sq  Mean Sq F value  Pr(>F)   
## SPP           1 0.00081 0.000812   0.827 0.36379   
## BASIN         6 0.02039 0.003398   3.463 0.00255 **
## SPP:BASIN     3 0.00256 0.000852   0.868 0.45782   
## Residuals   289 0.28360 0.000981                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ggplot(raw3, aes(x=factor(BASIN), y=abs.res.HL, fill=SPP)) +
  geom_boxplot()

A1.HD <- aov(abs.res.HD ~ SPP*BASIN, data=raw3)
summary(A1.HD)
##              Df  Sum Sq   Mean Sq F value Pr(>F)
## SPP           1 0.00000 0.0000000   0.000  0.999
## BASIN         6 0.00478 0.0007959   0.808  0.565
## SPP:BASIN     3 0.00382 0.0012743   1.293  0.277
## Residuals   289 0.28484 0.0009856
ggplot(raw3, aes(x=factor(BASIN), y=abs.res.HD, fill=SPP)) +
  geom_boxplot()

A1.HW <- aov(abs.res.HW ~ SPP*BASIN, data=raw3)
summary(A1.HW)
##              Df  Sum Sq   Mean Sq F value  Pr(>F)   
## SPP           1 0.00213 0.0021307   2.825 0.09387 . 
## BASIN         6 0.01505 0.0025085   3.326 0.00349 **
## SPP:BASIN     3 0.00076 0.0002544   0.337 0.79838   
## Residuals   289 0.21795 0.0007542                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ggplot(raw3, aes(x=factor(BASIN), y=abs.res.HW, fill=SPP)) +
  geom_boxplot()

A1.SnL <- aov(abs.res.SnL ~ SPP*BASIN, data=raw3)
summary(A1.SnL)
##              Df  Sum Sq  Mean Sq F value   Pr(>F)    
## SPP           1 0.00009 0.000095   0.162    0.688    
## BASIN         6 0.02054 0.003424   5.849 8.99e-06 ***
## SPP:BASIN     3 0.00238 0.000794   1.357    0.256    
## Residuals   289 0.16916 0.000585                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
A1.OL <- aov(abs.res.OL ~ SPP*BASIN, data=raw3)
summary(A1.OL)
##              Df    Sum Sq   Mean Sq F value Pr(>F)
## SPP           1 4.900e-33 4.915e-33   2.054  0.153
## BASIN         6 9.200e-33 1.535e-33   0.641  0.697
## SPP:BASIN     3 4.100e-33 1.358e-33   0.567  0.637
## Residuals   289 6.917e-31 2.394e-33
ggplot(raw3, aes(x=factor(BASIN), y=abs.res.OL, fill=SPP)) +
  geom_boxplot()

Watersheds
A2.D <- aov(abs.res.D ~ SPP*WATERSHED, data=raw3)
summary(A2.D)
##                Df Sum Sq Mean Sq F value Pr(>F)  
## SPP             1   0.06  0.0635   0.415 0.5198  
## WATERSHED      13   4.33  0.3331   2.178 0.0106 *
## SPP:WATERSHED   5   1.96  0.3922   2.564 0.0274 *
## Residuals     280  42.82  0.1529                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ggplot(raw3, aes(x=factor(WATERSHED), y=abs.res.D, fill=SPP)) +
  geom_boxplot()

A2.P1 <- aov(abs.res.P1 ~ SPP*WATERSHED, data=raw3)
summary(A2.P1)
##                Df Sum Sq Mean Sq F value Pr(>F)  
## SPP             1   0.29  0.2878   1.568 0.2116  
## WATERSHED      13   3.73  0.2868   1.562 0.0957 .
## SPP:WATERSHED   5   0.78  0.1569   0.855 0.5120  
## Residuals     280  51.41  0.1836                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
A2.P1.R <- aov(abs.res.P1.R ~ SPP*WATERSHED, data=raw3)
summary(A2.P1.R)
##                Df Sum Sq Mean Sq F value   Pr(>F)    
## SPP             1   1.80  1.7983  11.145 0.000957 ***
## WATERSHED      13   0.82  0.0629   0.390 0.972416    
## SPP:WATERSHED   5   0.16  0.0322   0.199 0.962521    
## Residuals     280  45.18  0.1614                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ggplot(raw3, aes(x=factor(WATERSHED), y=abs.res.P1.R, fill=SPP)) +
  geom_boxplot()

A2.LLSC <- aov(abs.res.LLSC ~ SPP*WATERSHED, data=raw3)
summary(A2.LLSC)
##                Df Sum Sq Mean Sq F value   Pr(>F)    
## SPP             1   1.40  1.4013   3.299   0.0704 .  
## WATERSHED      13  20.28  1.5596   3.672 2.15e-05 ***
## SPP:WATERSHED   5   4.75  0.9502   2.237   0.0509 .  
## Residuals     280 118.93  0.4248                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ggplot(raw3, aes(x=factor(WATERSHED), y=abs.res.LLSC, fill=SPP)) +
  geom_boxplot()

A2.SBLL <- aov(abs.res.SBLL ~ SPP*WATERSHED, data=raw3)
summary(A2.SBLL)
##                Df Sum Sq Mean Sq F value   Pr(>F)    
## SPP             1   0.65  0.6523   5.611   0.0185 *  
## WATERSHED      13   2.98  0.2289   1.969   0.0233 *  
## SPP:WATERSHED   5   4.57  0.9138   7.861 6.06e-07 ***
## Residuals     280  32.55  0.1162                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ggplot(raw3, aes(x=factor(WATERSHED), y=abs.res.SBLL, fill=SPP)) +
  geom_boxplot()

A2.BD <- aov(abs.res.BD ~ SPP*WATERSHED, data=raw3)
summary(A2.BD)
##                Df Sum Sq  Mean Sq F value Pr(>F)
## SPP             1 0.0023 0.002282   0.590  0.443
## WATERSHED      13 0.0310 0.002386   0.617  0.840
## SPP:WATERSHED   5 0.0275 0.005491   1.420  0.217
## Residuals     280 1.0829 0.003867
ggplot(raw3, aes(x=factor(WATERSHED), y=abs.res.BD, fill=SPP)) +
  geom_boxplot()

A2.CPD <- aov(abs.res.CPD ~ SPP*WATERSHED, data=raw3)
summary(A2.CPD)
##                Df  Sum Sq   Mean Sq F value Pr(>F)  
## SPP             1 0.00312 0.0031247   4.943 0.0270 *
## WATERSHED      13 0.01690 0.0012999   2.056 0.0169 *
## SPP:WATERSHED   5 0.00731 0.0014624   2.313 0.0441 *
## Residuals     280 0.17700 0.0006322                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ggplot(raw3, aes(x=factor(WATERSHED), y=abs.res.CPD, fill=SPP)) +
  geom_boxplot()

A2.CPL <- aov(abs.res.CPL ~ SPP*WATERSHED, data=raw3)
summary(A2.CPL)
##                Df  Sum Sq  Mean Sq F value   Pr(>F)    
## SPP             1 0.00000 0.000000   0.000   0.9866    
## WATERSHED      13 0.04566 0.003512   4.406 8.64e-07 ***
## SPP:WATERSHED   5 0.01021 0.002043   2.563   0.0275 *  
## Residuals     280 0.22322 0.000797                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
A2.PreDL <- aov(abs.res.PreDL ~ SPP*WATERSHED, data=raw3)
summary(A2.PreDL)
##                Df Sum Sq  Mean Sq F value Pr(>F)
## SPP             1 0.0012 0.001166   0.498  0.481
## WATERSHED      13 0.0220 0.001688   0.722  0.742
## SPP:WATERSHED   5 0.0128 0.002569   1.098  0.362
## Residuals     280 0.6552 0.002340
A2.DbL <- aov(abs.res.DbL ~ SPP*WATERSHED, data=raw3)
summary(A2.DbL)
##                Df  Sum Sq  Mean Sq F value   Pr(>F)    
## SPP             1 0.00556 0.005557   6.321   0.0125 *  
## WATERSHED      13 0.03949 0.003038   3.456 5.48e-05 ***
## SPP:WATERSHED   5 0.01333 0.002666   3.033   0.0110 *  
## Residuals     280 0.24615 0.000879                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ggplot(raw3, aes(x=factor(WATERSHED), y=abs.res.DbL, fill=SPP)) +
  geom_boxplot()

A2.HL <- aov(abs.res.HL ~ SPP*WATERSHED, data=raw3)
summary(A2.HL)
##                Df  Sum Sq  Mean Sq F value   Pr(>F)    
## SPP             1 0.00081 0.000812   0.896    0.345    
## WATERSHED      13 0.05036 0.003874   4.275 1.54e-06 ***
## SPP:WATERSHED   5 0.00247 0.000494   0.546    0.742    
## Residuals     280 0.25372 0.000906                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ggplot(raw3, aes(x=factor(WATERSHED), y=abs.res.HL, fill=SPP)) +
  geom_boxplot()

A2.HD <- aov(abs.res.HD ~ SPP*WATERSHED, data=raw3)
summary(A2.HD)
##                Df  Sum Sq   Mean Sq F value Pr(>F)
## SPP             1 0.00000 0.0000000   0.000  0.999
## WATERSHED      13 0.01725 0.0013271   1.371  0.173
## SPP:WATERSHED   5 0.00519 0.0010375   1.072  0.376
## Residuals     280 0.27100 0.0009678
ggplot(raw3, aes(x=factor(WATERSHED), y=abs.res.HD, fill=SPP)) +
  geom_boxplot()

A2.HW <- aov(abs.res.HW ~ SPP*WATERSHED, data=raw3)
summary(A2.HW)
##                Df  Sum Sq   Mean Sq F value Pr(>F)  
## SPP             1 0.00213 0.0021307   2.862 0.0918 .
## WATERSHED      13 0.02008 0.0015446   2.075 0.0157 *
## SPP:WATERSHED   5 0.00521 0.0010418   1.399 0.2246  
## Residuals     280 0.20848 0.0007446                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ggplot(raw3, aes(x=factor(WATERSHED), y=abs.res.HW, fill=SPP)) +
  geom_boxplot()

A2.SnL <- aov(abs.res.SnL ~ SPP*WATERSHED, data=raw3)
summary(A2.SnL)
##                Df  Sum Sq   Mean Sq F value   Pr(>F)    
## SPP             1 0.00009 0.0000948   0.180    0.672    
## WATERSHED      13 0.04027 0.0030976   5.890 1.23e-09 ***
## SPP:WATERSHED   5 0.00457 0.0009144   1.739    0.126    
## Residuals     280 0.14725 0.0005259                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
A2.OL <- aov(abs.res.OL ~ SPP*WATERSHED, data=raw3)
summary(A2.OL)
##                Df    Sum Sq   Mean Sq F value Pr(>F)
## SPP             1 4.900e-33 4.915e-33   2.037  0.155
## WATERSHED      13 2.900e-32 2.231e-33   0.925  0.528
## SPP:WATERSHED   5 2.000e-34 4.600e-35   0.019  1.000
## Residuals     280 6.758e-31 2.414e-33
ggplot(raw3, aes(x=factor(WATERSHED), y=abs.res.OL, fill=SPP)) +
  geom_boxplot()

Standardized

DON’T KNOW IF THIS IS USEFUL

The ANOVAs above focus on differences of particular traits as a factor of species and location. If we want to get an idea of variation in general as a factor of species and location, we can standardize the residuals (essentially unitless z-scores of residuals).

sd.res.D <- append(abs(sd.lat.D), abs(sd.form.D))
sd.res.P1 <- append(abs(sd.lat.P1), abs(sd.form.P1))
sd.res.P1.R <- append(abs(sd.lat.P1.R), abs(sd.form.P1.R))
sd.res.LLSC<- append(abs(sd.lat.LLSC), abs(sd.form.LLSC))
sd.res.SBLL<- append(abs(sd.lat.SBLL), abs(sd.form.SBLL))
sd.res.BD<- append(abs(sd.lat.BD), abs(sd.form.BD))
sd.res.CPD<- append(abs(sd.lat.CPD), abs(sd.form.CPD))
sd.res.CPL<- append(abs(sd.lat.CPL), abs(sd.form.CPL))
sd.res.PreDL <- append(abs(sd.lat.PreDL), abs(sd.form.PreDL))
sd.res.DbL <- append(abs(sd.lat.DbL), abs(sd.form.DbL))
sd.res.HL<- append(abs(sd.lat.HL), abs(sd.form.HL))
sd.res.HD<- append(abs(sd.lat.HD), abs(sd.form.HD))
sd.res.HW <- append(abs(sd.lat.HW), abs(sd.form.HW))
sd.res.SnL <- append(abs(sd.lat.SnL), abs(sd.form.SnL))
sd.res.OL <- append(abs(sd.lat.OL), abs(sd.form.OL))


raw4 <- cbind(raw3, sd.res.D, sd.res.P1, sd.res.P1.R, sd.res.LLSC, sd.res.SBLL, sd.res.BD, sd.res.CPD, sd.res.CPL, sd.res.PreDL, sd.res.DbL, sd.res.HL, sd.res.HD, sd.res.HW, sd.res.SnL, sd.res.OL)

raw5 <- cbind(raw4[1:14], stack(raw4[53:68])) 
## Warning in data.frame(..., check.names = FALSE): row names were found from a
## short variable and have been discarded
lat.raw5 <- raw5[raw5$SPP == "p.latipinna",]

form.raw5 <- raw5[raw5$SPP == "p.formosa",]

######ZONES#####

A3.lat <- aov(values~QUARTILE, data=lat.raw5)
summary(A3.lat)
##               Df Sum Sq Mean Sq F value Pr(>F)
## QUARTILE       3    0.3  0.1037   0.369  0.775
## Residuals   2140  601.4  0.2810
A3.form <- aov(values~QUARTILE, data=form.raw5)
summary(A3.form)
##               Df Sum Sq Mean Sq F value Pr(>F)
## QUARTILE       3    1.4  0.4703   1.579  0.192
## Residuals   2652  789.9  0.2979
#between species

A3 <- aov(values~QUARTILE*SPP, data=raw5)
summary(A3)
##                Df Sum Sq Mean Sq F value Pr(>F)  
## QUARTILE        3    0.6  0.2047   0.705 0.5489  
## SPP             1    1.0  1.0000   3.444 0.0635 .
## QUARTILE:SPP    3    1.1  0.3643   1.255 0.2882  
## Residuals    4792 1391.3  0.2903                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ggplot(raw5, aes(x=factor(QUARTILE), y=values, fill=SPP)) +
  geom_boxplot()

######BASINS#####

A4.lat <- aov(values~BASIN, data=lat.raw5)
summary(A4.lat)
##               Df Sum Sq Mean Sq F value Pr(>F)
## BASIN          5    1.1  0.2133   0.759  0.579
## Residuals   2138  600.6  0.2809
A4.form <- aov(values~BASIN, data=form.raw5)
summary(A4.form)
##               Df Sum Sq Mean Sq F value Pr(>F)
## BASIN          4    0.9  0.2307   0.774  0.542
## Residuals   2651  790.4  0.2982
#between species

A4 <- aov(values~BASIN*SPP, data=raw5)
summary(A4)
##               Df Sum Sq Mean Sq F value Pr(>F)  
## BASIN          6    2.0  0.3251   1.119 0.3482  
## SPP            1    0.8  0.7979   2.747 0.0975 .
## BASIN:SPP      3    0.2  0.0752   0.259 0.8551  
## Residuals   4789 1391.0  0.2905                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ggplot(raw5, aes(x=factor(BASIN), y=values, fill=SPP)) +
  geom_boxplot()

#####WATERSHEDS#####
A5.lat <- aov(values~WATERSHED, data=lat.raw5)
summary(A5.lat)
##               Df Sum Sq Mean Sq F value Pr(>F)  
## WATERSHED     11    6.1  0.5514   1.974 0.0273 *
## Residuals   2132  595.6  0.2794                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
A5.form <- aov(values~WATERSHED, data=form.raw5)
summary(A5.form)
##               Df Sum Sq Mean Sq F value Pr(>F)  
## WATERSHED      7    5.1  0.7267   2.447 0.0168 *
## Residuals   2648  786.2  0.2969                 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#between species

A5 <- aov(values~WATERSHED*SPP, data=raw5)
summary(A5)
##                 Df Sum Sq Mean Sq F value  Pr(>F)   
## WATERSHED       13    6.5  0.4965   1.717 0.05087 . 
## SPP              1    0.9  0.9387   3.247 0.07161 . 
## WATERSHED:SPP    5    4.7  0.9487   3.282 0.00583 **
## Residuals     4780 1381.8  0.2891                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
ggplot(raw5, aes(x=factor(WATERSHED), y=values, fill=SPP)) +
  geom_boxplot()

PCA analysis

PCA

LOGAN: CHECK THAT EACH VARIABLES IS NEAR NORMALLY DISTRIBUTED. IF NOT, LOG TRANSFORM BEFORE PCA. ALSO CHECK THAT PCA CALCULATES Z SCORES AND PLOTS BASED ON THAT; IF NOT CONVERT TO Z SCORES THEN RUN PCA.

In this analysis, I will compare the principle components after centering and scaling the data. A PCA analysis will help us determine what aspects of morphology influence the variation in our data the most without worrying about differences in scales/measurements. Currently, data consists of 116 Sailfin and 186 Amazon.

Chart

(chart <- matrix(c("Variable chart:", "D: dorsal ray count", "P1: left pectoral ray count", "P2.R: right pelvic rays", "P2.L: left pelvic rays", "A: anal rays", "P1.R: right pectoral ray count", "LLSC: lateral line scale count", "SALL: scales above lateral line", "SBLL: scales below lateral line", "SBDF: scales before dorsal fin", "TL: total length", "SL: standard length", "BD: body depth", "CPD: caudal peduncle depth", "CPL: caudal peduncle length", "PreDL: pre-dorsal length", "DbL: dorsal base length", "HL/HW/HD: head length/width/depth", "SnL: snout length", "OL: ocular length"), ncol=1, byrow=TRUE))
##       [,1]                               
##  [1,] "Variable chart:"                  
##  [2,] "D: dorsal ray count"              
##  [3,] "P1: left pectoral ray count"      
##  [4,] "P2.R: right pelvic rays"          
##  [5,] "P2.L: left pelvic rays"           
##  [6,] "A: anal rays"                     
##  [7,] "P1.R: right pectoral ray count"   
##  [8,] "LLSC: lateral line scale count"   
##  [9,] "SALL: scales above lateral line"  
## [10,] "SBLL: scales below lateral line"  
## [11,] "SBDF: scales before dorsal fin"   
## [12,] "TL: total length"                 
## [13,] "SL: standard length"              
## [14,] "BD: body depth"                   
## [15,] "CPD: caudal peduncle depth"       
## [16,] "CPL: caudal peduncle length"      
## [17,] "PreDL: pre-dorsal length"         
## [18,] "DbL: dorsal base length"          
## [19,] "HL/HW/HD: head length/width/depth"
## [20,] "SnL: snout length"                
## [21,] "OL: ocular length"

Loadings

raw1a <- subset(raw1, select = -c(16:18) )

PCA <- prcomp(raw1a[, 15:32], center=TRUE, scale. = TRUE) #includes all but pelvic traits, since they are the same in all species
summary(PCA)
## Importance of components:
##                          PC1    PC2     PC3     PC4     PC5     PC6     PC7
## Standard deviation     2.961 1.6259 1.08642 1.03324 0.97942 0.93729 0.80773
## Proportion of Variance 0.487 0.1469 0.06557 0.05931 0.05329 0.04881 0.03625
## Cumulative Proportion  0.487 0.6339 0.69942 0.75873 0.81203 0.86083 0.89708
##                            PC8     PC9    PC10    PC11   PC12    PC13    PC14
## Standard deviation     0.77182 0.54803 0.48623 0.46998 0.3746 0.31618 0.27529
## Proportion of Variance 0.03309 0.01669 0.01313 0.01227 0.0078 0.00555 0.00421
## Cumulative Proportion  0.93017 0.94686 0.95999 0.97226 0.9801 0.98561 0.98982
##                           PC15    PC16    PC17    PC18
## Standard deviation     0.24862 0.23425 0.22644 0.12339
## Proportion of Variance 0.00343 0.00305 0.00285 0.00085
## Cumulative Proportion  0.99326 0.99631 0.99915 1.00000
(loadings <- PCA$rotation[, 1:5])
##                PC1         PC2          PC3           PC4         PC5
## D     -0.025108263 -0.55543398 -0.204065448 -0.0432146833  0.01322127
## A     -0.007657967  0.05864241 -0.516347575 -0.3997522212 -0.61191359
## P1.R  -0.017685363  0.42264377 -0.121615925 -0.1355579521 -0.19302219
## LLSC  -0.078654354 -0.00632404 -0.479204155  0.5271970363  0.29628190
## SALL  -0.036439305 -0.05721360  0.519500405  0.4340355978 -0.54808415
## SBLL  -0.048429591 -0.05003098 -0.357637604  0.5491399767 -0.37340026
## SBDF  -0.004651513  0.52548030 -0.087495119  0.1774820105  0.14584307
## SL    -0.329799007  0.08491057  0.005044672 -0.0293415223  0.02269537
## BD    -0.320842190 -0.08571029 -0.005338963  0.0176425136  0.02068217
## CPD   -0.326190468 -0.02557070  0.047138928 -0.0583694362 -0.01044821
## CPL   -0.315318535  0.11009355  0.006579877 -0.0094553119  0.04336306
## PreDL -0.297122285  0.23349788  0.050939933 -0.0182778666 -0.02232858
## DbL   -0.245253913 -0.37105730 -0.047255019 -0.0675018348  0.01623612
## HL    -0.289172502  0.02061032 -0.069719784 -0.0865238724  0.04786649
## HD    -0.321400173 -0.04634301  0.022508677 -0.0496468119  0.06233383
## HW    -0.321926126  0.04469920  0.092090639  0.0326037283 -0.03172203
## SnL   -0.214460008 -0.05014512  0.130882266  0.0017726331 -0.16379453
## OL    -0.294285922  0.02932691 -0.061230019 -0.0001186146  0.03931992
sorted.loadings.1 <- loadings[order(loadings[, 1]), 1]
myTitle <- "Loadings Plot for PC1" 
myXlab  <- "Variable Loadings"
dotchart(sorted.loadings.1, main=myTitle, xlab=myXlab, cex=1.5, col="red")

sorted.loadings.2 <- loadings[order(loadings[, 2]), 2]
myTitle <- "Loadings Plot for PC2" 
myXlab  <- "Variable Loadings"
dotchart(sorted.loadings.2, main=myTitle, xlab=myXlab, cex=1.5, col="red")

sorted.loadings.3 <- loadings[order(loadings[, 3]), 3]
myTitle <- "Loadings Plot for PC3" 
myXlab  <- "Variable Loadings"
dotchart(sorted.loadings.3, main=myTitle, xlab=myXlab, cex=1.5, col="red")

VM_PCA <- varimax(PCA$rotation) 

VM_PCA
## $loadings
## 
## Loadings:
##       PC1 PC2 PC3 PC4 PC5 PC6 PC7 PC8 PC9 PC10 PC11 PC12 PC13 PC14 PC15 PC16
## D         -1                                                                
## A                     -1                                                    
## P1.R                          -1                                            
## LLSC               1                                                        
## SALL           1                                                            
## SBLL                      -1                                                
## SBDF                                       1                                
## SL                                                                          
## BD                                                                      -1  
## CPD                                                                         
## CPL                                                 -1                      
## PreDL                                                              -1       
## DbL                                                       1                 
## HL                                     1                                    
## HD                                                             1            
## HW    -1                                                                    
## SnL                               -1                                        
## OL                                             -1                           
##       PC17 PC18
## D              
## A              
## P1.R           
## LLSC           
## SALL           
## SBLL           
## SBDF           
## SL          1  
## BD             
## CPD    1       
## CPL            
## PreDL          
## DbL            
## HL             
## HD             
## HW             
## SnL            
## OL             
## 
##                  PC1   PC2   PC3   PC4   PC5   PC6   PC7   PC8   PC9  PC10
## SS loadings    1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000
## Proportion Var 0.056 0.056 0.056 0.056 0.056 0.056 0.056 0.056 0.056 0.056
## Cumulative Var 0.056 0.111 0.167 0.222 0.278 0.333 0.389 0.444 0.500 0.556
##                 PC11  PC12  PC13  PC14  PC15  PC16  PC17  PC18
## SS loadings    1.000 1.000 1.000 1.000 1.000 1.000 1.000 1.000
## Proportion Var 0.056 0.056 0.056 0.056 0.056 0.056 0.056 0.056
## Cumulative Var 0.611 0.667 0.722 0.778 0.833 0.889 0.944 1.000
## 
## $rotmat
##                [,1]        [,2]         [,3]         [,4]         [,5]
##  [1,]  0.3219261237  0.02510828 -0.036439306 -0.078654368  0.007657989
##  [2,] -0.0446992069  0.55543394 -0.057213641 -0.006324052 -0.058642377
##  [3,] -0.0920906421  0.20406545  0.519500670 -0.479204308  0.516347268
##  [4,] -0.0326037378  0.04321468  0.434035867  0.527197159  0.399751969
##  [5,]  0.0317220321 -0.01322128 -0.548083772  0.296281819  0.611913921
##  [6,]  0.0136429476 -0.04394442  0.388233014  0.529602504 -0.351562759
##  [7,]  0.0300222480  0.15540493 -0.155600536 -0.172681388 -0.193112792
##  [8,] -0.0263614684  0.04528755  0.164534883 -0.167951397 -0.124539401
##  [9,]  0.1342223346 -0.08204914  0.087627176  0.012441859  0.034150254
## [10,] -0.1036782848 -0.44030311  0.074088336 -0.172666056  0.006036049
## [11,] -0.2364987676  0.26371383 -0.126067295  0.092162379 -0.115021117
## [12,] -0.0465619994  0.53501853 -0.002276466  0.093586039 -0.048447786
## [13,]  0.4763014542  0.21639007  0.053353969 -0.039597216 -0.018100105
## [14,]  0.0004626634 -0.04444238  0.026277094  0.010977539 -0.016082025
## [15,] -0.0928579707 -0.05365086  0.003045335 -0.098490444 -0.034010892
## [16,] -0.5107289452 -0.05121511 -0.025888687  0.027609857 -0.024839336
## [17,]  0.5343922581 -0.09062037  0.033283349  0.023118943  0.003021833
## [18,]  0.1151678784  0.01195750  0.018346326  0.002323942  0.009991313
##               [,6]        [,7]         [,8]        [,9]         [,10]
##  [1,]  0.048429569  0.01768536  0.214460009 -0.28917250 -0.0046515118
##  [2,]  0.050030984 -0.42264377  0.050145118  0.02061032  0.5254803108
##  [3,]  0.357637452  0.12161593 -0.130882266 -0.06971978 -0.0874951163
##  [4,] -0.549139829  0.13555795 -0.001772633 -0.08652387  0.1774820114
##  [5,]  0.373400344  0.19302219  0.163794533  0.04786649  0.1458430706
##  [6,]  0.635932024  0.15025751 -0.038190321 -0.02887285  0.0738524609
##  [7,] -0.110709157  0.74061288 -0.376784024  0.02128196  0.3299828233
##  [8,] -0.051549852  0.31444972  0.866326081  0.15705131  0.1017346315
##  [9,] -0.005265660 -0.01312401 -0.088985564  0.80807722  0.1093110673
## [10,]  0.060467289 -0.13242476 -0.002448353 -0.22909667  0.6987442332
## [11,] -0.035905583  0.20660581  0.055806882 -0.30184685 -0.1331149106
## [12,]  0.017701395 -0.08227240 -0.007200975  0.05666880 -0.0053772537
## [13,] -0.016079039  0.08390788 -0.021282779  0.00476731  0.1277171583
## [14,] -0.023990107 -0.02159300 -0.007990821 -0.22466796  0.0182482439
## [15,] -0.020446257  0.03077480 -0.017531406  0.02720849  0.0477060736
## [16,]  0.001789952  0.07105234  0.040017163  0.07723549  0.0329101541
## [17,] -0.019811963  0.00238789  0.006509850 -0.13359655 -0.0008947563
## [18,] -0.002592149 -0.00092466  0.005121886 -0.02270106 -0.0280220293
##               [,11]         [,12]       [,13]       [,14]         [,15]
##  [1,]  0.2942859234  0.3153192664 -0.24525297 -0.32140019  0.2971222926
##  [2,] -0.0293269102 -0.1100924591 -0.37105766 -0.04634301 -0.2334978927
##  [3,]  0.0612300201 -0.0065797470 -0.04725504  0.02250868 -0.0509399352
##  [4,]  0.0001186154  0.0094555097 -0.06750181 -0.04964681  0.0182778668
##  [5,] -0.0393199164 -0.0433631099  0.01623598  0.06233384  0.0223285794
##  [6,] -0.0324815176 -0.0190456592 -0.05686694 -0.02430586  0.0006229111
##  [7,] -0.1366527484  0.1489676508 -0.14611544  0.03068266 -0.0646743205
##  [8,] -0.1130938511 -0.0004479098  0.01634738  0.11054770 -0.0852014182
##  [9,]  0.5008070786  0.0580593824 -0.01493268  0.11738705  0.0245567369
## [10,]  0.1981464875  0.1331502612  0.30269250  0.03103419  0.1110463216
## [11,]  0.7557661500 -0.1778142641  0.07766334  0.02973031 -0.1363150999
## [12,] -0.0696021143  0.3843471014  0.36505820  0.29552124  0.4728614204
## [13,] -0.0341062726 -0.5167317622  0.52407171 -0.30400581  0.1226784030
## [14,]  0.0109614331 -0.1210235434  0.10846697  0.70825624 -0.0355987528
## [15,] -0.0070553388 -0.5203739096 -0.35389964  0.16845443  0.6281445377
## [16,] -0.0242978744 -0.1076411066 -0.07390241 -0.24068768  0.3092391095
## [17,]  0.0772508602 -0.0390618306 -0.33669606  0.29203155  0.0225063642
## [18,]  0.0270806510  0.3128669202 -0.09290549 -0.07313309  0.2831870115
##              [,16]        [,17]        [,18]
##  [1,]  0.320842156 -0.326190480 -0.329799006
##  [2,]  0.085710295 -0.025570694  0.084910572
##  [3,]  0.005338974  0.047138935  0.005044672
##  [4,] -0.017642508 -0.058369435 -0.029341522
##  [5,] -0.020682173 -0.010448216  0.022695375
##  [6,] -0.010450947 -0.054650898 -0.003282193
##  [7,]  0.035665976 -0.046116425 -0.028552797
##  [8,] -0.063816178  0.022817476  0.047558720
##  [9,]  0.135533534 -0.069632230 -0.034558486
## [10,] -0.167139303  0.090000176 -0.013181473
## [11,] -0.169309096  0.099637999  0.118722245
## [12,] -0.155683617  0.152330270 -0.214231096
## [13,]  0.186265919  0.034812983  0.106212379
## [14,]  0.579174916 -0.266359247  0.096040493
## [15,] -0.311068254 -0.238586348 -0.015193539
## [16,]  0.561394466  0.483048937 -0.032924944
## [17,] -0.035577532  0.687773168 -0.072990532
## [18,]  0.030876634 -0.006467397  0.889397637

Plots

Component extractions
library(AMR) 
library(ggplot2)
library(ggfortify)

(PCA$sdev ^ 2)
##  [1] 8.76586176 2.64347463 1.18030344 1.06757490 0.95926986 0.87851453
##  [7] 0.65242074 0.59570413 0.30033241 0.23641483 0.22088004 0.14031638
## [13] 0.09996690 0.07578261 0.06181091 0.05487148 0.05127637 0.01522408
evplot <- function(ev)
{
    # Broken stick model (MacArthur 1957)
    n <- length(ev)
    bsm <- data.frame(j=seq(1:n), p=0)
    bsm$p[1] <- 1/n
    for (i in 2:n) bsm$p[i] <- bsm$p[i-1] + (1/(n + 1 - i))
    bsm$p <- 100*bsm$p/n
    # Plot eigenvalues and % of variation for each axis
    op <- par(mfrow=c(2,1))
    barplot(ev, main="Eigenvalues", col="bisque", las=2)
    abline(h=mean(ev), col="red")
    legend("topright", "Average eigenvalue", lwd=1, col=2, bty="n")
    barplot(t(cbind(100*ev/sum(ev), bsm$p[n:1])), beside=TRUE, 
        main="% variation", col=c("bisque",2), las=2)
    legend("topright", c("% eigenvalue", "Broken stick model"), 
        pch=15, col=c("bisque",2), bty="n")
    par(op)
}


ev <- PCA$sdev^2 
evplot(ev) #according to Kaiser-Guttman criteron, we can use the first 4 PCs, even though the broken stick model shows only the first above the red bar plot... not 100% confident I know what this means, but pretty sure PC1 is body size

PCA 1v2
plot1<- autoplot(PCA, data = raw1, colour='SPP', loadings=FALSE, loadings.label=FALSE, frame=TRUE, frame.type='norm')+ ggtitle("PCA Plot of Morphology traits") + theme_minimal() 
plot1

By zone
plot2<- autoplot(PCA, data = raw1, colour='QUARTILE', shape="SPP", frame=TRUE, frame.type='norm')+ ggtitle("PCA Plot of Morphology traits") + theme_minimal() 
plot2

By Basin
plot3<- autoplot(PCA, data = raw1, colour='BASIN', shape="SPP", frame=TRUE, frame.type='norm')+ ggtitle("PCA Plot of Morphology traits") + theme_minimal() 
plot3

By watershed
plot4<- autoplot(PCA, data = raw1, colour='WATERSHED', shape="SPP", frame=TRUE, frame.type='norm')+ ggtitle("PCA Plot of Morphology traits") + theme_minimal() 
plot4
## Too few points to calculate an ellipse
## Too few points to calculate an ellipse

PCA 2v3
plot5<- autoplot(PCA, x=2, y=3, data = raw1, colour='SPP', loadings=FALSE, loadings.label=FALSE, frame=TRUE, frame.type='norm')+ ggtitle("PCA Plot of Morphology traits") + theme_minimal()
plot5

By zone
plot6<- autoplot(PCA, x=2, y=3, data = raw1, colour='QUARTILE', shape="SPP", frame=TRUE, frame.type='norm')+ ggtitle("PCA Plot of Morphology traits") + theme_minimal() 
plot6

By basin
plot7<- autoplot(PCA, x=2, y=3, data = raw1, colour='BASIN', shape="SPP", frame=TRUE, frame.type='norm')+ ggtitle("PCA Plot of Morphology traits") + theme_minimal() 
plot7

By watershed
plot8<- autoplot(PCA, x=2, y=3, data = raw1, colour='WATERSHED', shape="SPP", frame=TRUE, frame.type='norm')+ ggtitle("PCA Plot of Morphology traits") + theme_minimal() 
plot8
## Too few points to calculate an ellipse
## Too few points to calculate an ellipse







Filtered data

Now I will repeat analysis on a dataset that filters for fish above 32.5mm (this is the average for the range I found in literature, 28-37mm). This would confirm that all fish are likely adults, and remove any bias from juvenile proportions (ex. if juvenile head is proportionately larger than expected for body).

Initial steps

Data collection

library(dplyr)

rawF <- filter(raw1,SL>=32.5)

Checking normality

Checking normality

I will use Shapiro-wilke, histograms, and QQ plots to determine what traits are normal. These will only be performed on continuous variables, as discrete variables are not normal by nature.

Conclusions: literally all of them are NOT normal… will log transform them and run parametric tests (t-test, F-test, ANOVA). Not sure what to do with the discrete variables…[Logan said to use non-parametric OR a glmer with poisson distribution]

SW Test
shapiro.test(rawF$SL)
## 
##  Shapiro-Wilk normality test
## 
## data:  rawF$SL
## W = 0.88785, p-value = 9.323e-12
shapiro.test(rawF$BD)
## 
##  Shapiro-Wilk normality test
## 
## data:  rawF$BD
## W = 0.94379, p-value = 1.537e-07
shapiro.test(rawF$CPD)
## 
##  Shapiro-Wilk normality test
## 
## data:  rawF$CPD
## W = 0.9366, p-value = 3.381e-08
shapiro.test(rawF$CPL)
## 
##  Shapiro-Wilk normality test
## 
## data:  rawF$CPL
## W = 0.91011, p-value = 2.741e-10
shapiro.test(rawF$PreDL)
## 
##  Shapiro-Wilk normality test
## 
## data:  rawF$PreDL
## W = 0.94232, p-value = 1.118e-07
shapiro.test(rawF$DbL)
## 
##  Shapiro-Wilk normality test
## 
## data:  rawF$DbL
## W = 0.98676, p-value = 0.03793
shapiro.test(rawF$HL)
## 
##  Shapiro-Wilk normality test
## 
## data:  rawF$HL
## W = 0.87215, p-value = 1.115e-12
shapiro.test(rawF$HD)
## 
##  Shapiro-Wilk normality test
## 
## data:  rawF$HD
## W = 0.91737, p-value = 9.289e-10
shapiro.test(rawF$HW)
## 
##  Shapiro-Wilk normality test
## 
## data:  rawF$HW
## W = 0.92699, p-value = 5.211e-09
shapiro.test(rawF$SnL)
## 
##  Shapiro-Wilk normality test
## 
## data:  rawF$SnL
## W = 0.56341, p-value < 2.2e-16
shapiro.test(rawF$OL)
## 
##  Shapiro-Wilk normality test
## 
## data:  rawF$OL
## W = 0.9821, p-value = 0.006707
Histograms
hist(rawF$SL)

hist(rawF$BD)

hist(rawF$CPD)

hist(rawF$CPL)

hist(rawF$PreDL)

hist(rawF$DbL)

hist(rawF$HL)

hist(rawF$HD)

hist(rawF$HW)

hist(rawF$SnL)

hist(rawF$OL)

QQ plots
qqnorm(rawF$SL)
qqline(rawF$SL)

qqnorm(rawF$BD)
qqline(rawF$BD)

qqnorm(rawF$CPD)
qqline(rawF$CPD)

qqnorm(rawF$CPL)
qqline(rawF$CPL)

qqnorm(rawF$PreDL)
qqline(rawF$PreDL)

qqnorm(rawF$DbL)
qqline(rawF$DbL)

qqnorm(rawF$HL)
qqline(rawF$HL)

qqnorm(rawF$HD)
qqline(rawF$HD)

qqnorm(rawF$HW)
qqline(rawF$HW)

qqnorm(rawF$SnL)
qqline(rawF$SnL)

qqnorm(rawF$OL)
qqline(rawF$OL)

Log transformations

Since all of the continuous characters were not normal, I will log transform them, and proceed with the transformed data in all analyses.

rawF[paste0(names(rawF), '_log')] <- log(rawF[, 25:35], 10)
#for some reason did the log of the whole data frame instead of those specific columns, so I will just remove the unnecessary columns instead of spending time on what went wrong.

rawF <- subset(rawF, select = -c(48:77) )

#double checking normality with a SW test and QQ plot

shapiro.test(rawF$SL_log)
## 
##  Shapiro-Wilk normality test
## 
## data:  rawF$SL_log
## W = 0.98318, p-value = 0.009921
shapiro.test(rawF$BD_log)
## 
##  Shapiro-Wilk normality test
## 
## data:  rawF$BD_log
## W = 0.95211, p-value = 1.024e-06
shapiro.test(rawF$CPD_log)
## 
##  Shapiro-Wilk normality test
## 
## data:  rawF$CPD_log
## W = 0.97428, p-value = 0.000465
shapiro.test(rawF$CPL_log)
## 
##  Shapiro-Wilk normality test
## 
## data:  rawF$CPL_log
## W = 0.97841, p-value = 0.001827
shapiro.test(rawF$PreDL_log)
## 
##  Shapiro-Wilk normality test
## 
## data:  rawF$PreDL_log
## W = 0.63075, p-value < 2.2e-16
shapiro.test(rawF$DbL_log)
## 
##  Shapiro-Wilk normality test
## 
## data:  rawF$DbL_log
## W = 0.97229, p-value = 0.0002479
shapiro.test(rawF$HL_log)
## 
##  Shapiro-Wilk normality test
## 
## data:  rawF$HL_log
## W = 0.97621, p-value = 0.0008721
shapiro.test(rawF$HD_log)
## 
##  Shapiro-Wilk normality test
## 
## data:  rawF$HD_log
## W = 0.88455, p-value = 5.875e-12
shapiro.test(rawF$HW_log)
## 
##  Shapiro-Wilk normality test
## 
## data:  rawF$HW_log
## W = 0.98022, p-value = 0.003423
shapiro.test(rawF$SnL_log)
## 
##  Shapiro-Wilk normality test
## 
## data:  rawF$SnL_log
## W = 0.93195, p-value = 1.341e-08
shapiro.test(rawF$OL_log)
## 
##  Shapiro-Wilk normality test
## 
## data:  rawF$OL_log
## W = 0.98727, p-value = 0.04594
qqnorm(rawF$SL_log)
qqline(rawF$SL_log)

qqnorm(rawF$BD_log)
qqline(rawF$BD_log)

qqnorm(rawF$CPD_log)
qqline(rawF$CPD_log)

qqnorm(rawF$CPL_log)
qqline(rawF$CPL_log)

qqnorm(rawF$PreDL_log)
qqline(rawF$PreDL_log)

qqnorm(rawF$DbL_log)
qqline(rawF$DbL_log)

qqnorm(rawF$HL_log)
qqline(rawF$HL_log)

qqnorm(rawF$HD_log)
qqline(rawF$HD_log)

qqnorm(rawF$HW_log)
qqline(rawF$HW_log)

qqnorm(rawF$SnL_log)
qqline(rawF$SnL_log)

qqnorm(rawF$OL_log)
qqline(rawF$OL_log)

Correcting for body size (residuals)

SL correction

Since amazons are in general bigger than sailfin, we don’t want any results to be due to this difference in body size bias. Therefore, we will see what traits are influenced by body size (regressions) and correct for body size when necessary (absolute value of residuals). We can then use the residuals when comparing between species for traits that are influenced by body size, and raw/transformed data for traits that are not influenced by body size. I will also calculate standardized residuals to compare residuals across traits in later analyses.

Quick results summary: traits not influenced by body size are left & right pelvic, anal, scales above and below lateral line (except in mexicana), scales before dorsal fin and fluctuating asymmetry; all other traits influenced by body size.

Regressions

  • STEP ONE: plot trait vs body size
library(ggplot2)
library(ggpubr)



lat.F <- rawF[rawF$SPP == "p.latipinna",]


form.F <- rawF[rawF$SPP == "p.formosa",]


mex.F <- rawF[rawF$SPP == "p.mexicana",]



##### LAT #####
F.reg.lat.D <- lm(lat.F$D ~ lat.F$SL)
F.sd.lat.D <- rstandard(F.reg.lat.D)
F.reg.lat.D.plot <- ggplot(lat.F, aes(x = SL, y = D)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
F.reg.lat.D.plot
## `geom_smooth()` using formula = 'y ~ x'

F.reg.lat.P1 <- lm(lat.F$P1 ~ lat.F$SL)
F.sd.lat.P1 <- rstandard(F.reg.lat.P1)
F.reg.lat.P1.plot <- ggplot(lat.F, aes(x = SL, y = P1)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
F.reg.lat.P1.plot
## `geom_smooth()` using formula = 'y ~ x'

F.reg.lat.P2.L <- lm(lat.F$P2.L ~ lat.F$SL)
F.sd.lat.P2.L <- rstandard(F.reg.lat.P2.L)
F.reg.lat.P2.L.plot <- ggplot(lat.F, aes(x = SL, y = P2.L)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
F.reg.lat.P2.L.plot
## `geom_smooth()` using formula = 'y ~ x'

F.reg.lat.P2.R <- lm(lat.F$P2.R ~ lat.F$SL)
F.sd.lat.P2.R <- rstandard(F.reg.lat.P2.R)
F.reg.lat.P2.R.plot <- ggplot(lat.F, aes(x = SL, y = P2.R)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
F.reg.lat.P2.R.plot
## `geom_smooth()` using formula = 'y ~ x'

F.reg.lat.A <- lm(lat.F$A ~ lat.F$SL)
F.sd.lat.A <- rstandard(F.reg.lat.A)
F.reg.lat.A.plot <- ggplot(lat.F, aes(x = SL, y = A)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
F.reg.lat.A.plot
## `geom_smooth()` using formula = 'y ~ x'

F.reg.lat.P1.R <- lm(lat.F$P1.R ~ lat.F$SL)
F.sd.lat.P1.R <- rstandard(F.reg.lat.P1.R)
F.reg.lat.P1.R.plot <- ggplot(lat.F, aes(x = SL, y = P1.R)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
F.reg.lat.P1.R.plot
## `geom_smooth()` using formula = 'y ~ x'

F.reg.lat.LLSC <- lm(lat.F$LLSC ~ lat.F$SL)
F.sd.lat.LLSC <- rstandard(F.reg.lat.LLSC)
F.reg.lat.LLSC.plot <- ggplot(lat.F, aes(x = SL, y = LLSC)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
F.reg.lat.LLSC.plot
## `geom_smooth()` using formula = 'y ~ x'

F.reg.lat.SALL <- lm(lat.F$SALL ~ lat.F$SL)
F.sd.lat.SALL <- rstandard(F.reg.lat.SALL)
F.reg.lat.SALL.plot <- ggplot(lat.F, aes(x = SL, y = SALL)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
F.reg.lat.SALL.plot
## `geom_smooth()` using formula = 'y ~ x'

F.reg.lat.SBLL <- lm(lat.F$SBLL ~ lat.F$SL)
F.sd.lat.SBLL <- rstandard(F.reg.lat.SBLL)
F.reg.lat.SBLL.plot <- ggplot(lat.F, aes(x = SL, y = SBLL)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
F.reg.lat.SBLL.plot
## `geom_smooth()` using formula = 'y ~ x'

F.reg.lat.SBDF <- lm(lat.F$SBDF ~ lat.F$SL)
F.sd.lat.SBDF <- rstandard(F.reg.lat.SBDF)
F.reg.lat.SBDF.plot <- ggplot(lat.F, aes(x = SL, y = SBDF)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
F.reg.lat.SBDF.plot
## `geom_smooth()` using formula = 'y ~ x'

F.reg.lat.BD <- lm(lat.F$BD_log ~ lat.F$SL_log)
F.sd.lat.BD <- rstandard(F.reg.lat.BD)
F.reg.lat.BD.plot <- ggplot(lat.F, aes(x = SL_log, y = BD_log)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
F.reg.lat.BD.plot
## `geom_smooth()` using formula = 'y ~ x'

F.reg.lat.CPD <- lm(lat.F$CPD_log ~ lat.F$SL_log)
F.sd.lat.CPD <- rstandard(F.reg.lat.CPD)
F.reg.lat.CPD.plot <- ggplot(lat.F, aes(x = SL_log, y = CPD_log)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
F.reg.lat.CPD.plot
## `geom_smooth()` using formula = 'y ~ x'

F.reg.lat.CPL <- lm(lat.F$CPL_log ~ lat.F$SL_log)
F.sd.lat.CPL <- rstandard(F.reg.lat.CPL)
F.reg.lat.CPL.plot <- ggplot(lat.F, aes(x = SL_log, y = CPL_log)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
F.reg.lat.CPL.plot
## `geom_smooth()` using formula = 'y ~ x'

F.reg.lat.PreDL <- lm(lat.F$PreDL_log ~ lat.F$SL_log)
F.sd.lat.PreDL <- rstandard(F.reg.lat.PreDL)
F.reg.lat.PreDL.plot <- ggplot(lat.F, aes(x = SL_log, y = PreDL_log)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
F.reg.lat.PreDL.plot
## `geom_smooth()` using formula = 'y ~ x'

F.reg.lat.DbL <- lm(lat.F$DbL_log ~ lat.F$SL_log)
F.sd.lat.DbL <- rstandard(F.reg.lat.DbL)
F.reg.lat.DbL.plot <- ggplot(lat.F, aes(x = SL_log, y = DbL_log)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
F.reg.lat.DbL.plot
## `geom_smooth()` using formula = 'y ~ x'

F.reg.lat.HL <- lm(lat.F$HL_log ~ lat.F$SL_log)
F.sd.lat.HL <- rstandard(F.reg.lat.HL)
F.reg.lat.HL.plot <- ggplot(lat.F, aes(x = SL_log, y = HL_log)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
F.reg.lat.HL.plot
## `geom_smooth()` using formula = 'y ~ x'

F.reg.lat.HD <- lm(lat.F$HD_log ~ lat.F$SL_log)
F.sd.lat.HD <- rstandard(F.reg.lat.HD)
F.reg.lat.HD.plot <- ggplot(lat.F, aes(x = SL_log, y = HD_log)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
F.reg.lat.HD.plot
## `geom_smooth()` using formula = 'y ~ x'

F.reg.lat.HW <- lm(lat.F$HW_log ~ lat.F$SL_log)
F.sd.lat.HW <- rstandard(F.reg.lat.HW)
F.reg.lat.HW.plot <- ggplot(lat.F, aes(x = SL_log, y = HW_log)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
F.reg.lat.HW.plot
## `geom_smooth()` using formula = 'y ~ x'

F.reg.lat.SnL <- lm(lat.F$SnL_log ~ lat.F$SL_log)
F.sd.lat.SnL <- rstandard(F.reg.lat.SnL)
F.reg.lat.SnL.plot <- ggplot(lat.F, aes(x = SL_log, y = SnL_log)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
F.reg.lat.SnL.plot
## `geom_smooth()` using formula = 'y ~ x'

F.reg.lat.OL <- lm(lat.F$OL_log ~ lat.F$SL_log)
F.sd.lat.OL <- rstandard(F.reg.lat.OL)
F.reg.lat.OL.plot <- ggplot(lat.F, aes(x = SL_log, y = OL_log)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
F.reg.lat.OL.plot
## `geom_smooth()` using formula = 'y ~ x'

F.reg.lat.FLA <- lm(lat.F$FLA ~ lat.F$SL)
F.sd.lat.FLA <- rstandard(F.reg.lat.FLA)
F.reg.lat.FLA.plot <- ggplot(lat.F, aes(x = SL, y = FLA)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
F.reg.lat.FLA.plot
## `geom_smooth()` using formula = 'y ~ x'

##### FORM #####

F.reg.form.D <- lm(form.F$D ~ form.F$SL)
F.sd.form.D <- rstandard(F.reg.form.D)
F.reg.form.D.plot <- ggplot(form.F, aes(x =SL, y = D)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
F.reg.form.D.plot
## `geom_smooth()` using formula = 'y ~ x'

F.reg.form.P1 <- lm(form.F$P1 ~ form.F$SL)
F.sd.form.P1 <- rstandard(F.reg.form.P1)
F.reg.form.P1.plot <- ggplot(form.F, aes(x = SL, y = P1)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
F.reg.form.P1.plot
## `geom_smooth()` using formula = 'y ~ x'

F.reg.form.P2.L <- lm(form.F$P2.L ~ form.F$SL)
F.sd.form.P2.L <- rstandard(F.reg.form.P2.L)
F.reg.form.P2.L.plot <- ggplot(form.F, aes(x = SL, y = P2.L)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
F.reg.form.P2.L.plot
## `geom_smooth()` using formula = 'y ~ x'

F.reg.form.P2.R <- lm(form.F$P2.R ~ form.F$SL)
F.sd.form.P2.R <- rstandard(F.reg.form.P2.R)
F.reg.form.P2.R.plot <- ggplot(form.F, aes(x = SL, y = P2.R)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
F.reg.form.P2.R.plot
## `geom_smooth()` using formula = 'y ~ x'

F.reg.form.A <- lm(form.F$A ~ form.F$SL)
F.sd.form.A <- rstandard(F.reg.form.A)
F.reg.form.A.plot <- ggplot(form.F, aes(x = SL, y = A)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
F.reg.form.A.plot
## `geom_smooth()` using formula = 'y ~ x'

F.reg.form.P1.R <- lm(form.F$P1.R ~ form.F$SL)
F.sd.form.P1.R <- rstandard(F.reg.form.P1.R)
F.reg.form.P1.R.plot <- ggplot(form.F, aes(x = SL, y = P1.R)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
F.reg.form.P1.R.plot
## `geom_smooth()` using formula = 'y ~ x'

F.reg.form.LLSC <- lm(form.F$LLSC ~ form.F$SL)
F.sd.form.LLSC <- rstandard(F.reg.form.LLSC)
F.reg.form.LLSC.plot <- ggplot(form.F, aes(x = SL, y = LLSC)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
F.reg.form.LLSC.plot
## `geom_smooth()` using formula = 'y ~ x'

F.reg.form.SALL <- lm(form.F$SALL ~ form.F$SL)
F.sd.form.SALL <- rstandard(F.reg.form.SALL)
F.reg.form.SALL.plot <- ggplot(form.F, aes(x = SL, y = SALL)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
F.reg.form.SALL.plot
## `geom_smooth()` using formula = 'y ~ x'

F.reg.form.SBLL <- lm(form.F$SBLL ~ form.F$SL)
F.sd.form.SBLL <- rstandard(F.reg.form.SBLL)
F.reg.form.SBLL.plot <- ggplot(form.F, aes(x = SL, y = SBLL)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
F.reg.form.SBLL.plot
## `geom_smooth()` using formula = 'y ~ x'

F.reg.form.SBDF <- lm(form.F$SBDF ~ form.F$SL)
F.sd.form.SBDF <- rstandard(F.reg.form.SBDF)
F.reg.form.SBDF.plot <- ggplot(form.F, aes(x = SL, y = SBDF)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
F.reg.form.SBDF.plot
## `geom_smooth()` using formula = 'y ~ x'

F.reg.form.BD <- lm(form.F$BD_log ~ form.F$SL_log)
F.sd.form.BD <- rstandard(F.reg.form.BD)
F.reg.form.BD.plot <- ggplot(form.F, aes(x = SL_log, y = BD_log)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
F.reg.form.BD.plot
## `geom_smooth()` using formula = 'y ~ x'

F.reg.form.CPD <- lm(form.F$CPD_log ~ form.F$SL_log)
F.sd.form.CPD <- rstandard(F.reg.form.CPD)
F.reg.form.CPD.plot <- ggplot(form.F, aes(x = SL_log, y = CPD_log)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
F.reg.form.CPD.plot
## `geom_smooth()` using formula = 'y ~ x'

F.reg.form.CPL <- lm(form.F$CPL_log ~ form.F$SL_log)
F.sd.form.CPL <- rstandard(F.reg.form.CPL)
F.reg.form.CPL.plot <- ggplot(form.F, aes(x = SL_log, y = CPL_log)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
F.reg.form.CPL.plot
## `geom_smooth()` using formula = 'y ~ x'

F.reg.form.PreDL <- lm(form.F$PreDL_log ~ form.F$SL_log)
F.sd.form.PreDL <- rstandard(F.reg.form.PreDL)
F.reg.form.PreDL.plot <- ggplot(form.F, aes(x = SL_log, y = PreDL_log)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
F.reg.form.PreDL.plot
## `geom_smooth()` using formula = 'y ~ x'

F.reg.form.DbL <- lm(form.F$DbL_log ~ form.F$SL_log)
F.sd.form.DbL <- rstandard(F.reg.form.DbL)
F.reg.form.DbL.plot <- ggplot(form.F, aes(x = SL_log, y = DbL_log)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
F.reg.form.DbL.plot
## `geom_smooth()` using formula = 'y ~ x'

F.reg.form.HL <- lm(form.F$HL_log ~ form.F$SL_log)
F.sd.form.HL <- rstandard(F.reg.form.HL)
F.reg.form.HL.plot <- ggplot(form.F, aes(x = SL_log, y = HL_log)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
F.reg.form.HL.plot
## `geom_smooth()` using formula = 'y ~ x'

F.reg.form.HD <- lm(form.F$HD_log ~ form.F$SL_log)
F.sd.form.HD <- rstandard(F.reg.form.HD)
F.reg.form.HD.plot <- ggplot(form.F, aes(x = SL_log, y = HD_log)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
F.reg.form.HD.plot
## `geom_smooth()` using formula = 'y ~ x'

F.reg.form.HW <- lm(form.F$HW_log ~ form.F$SL_log)
F.sd.form.HW <- rstandard(F.reg.form.HW)
F.reg.form.HW.plot <- ggplot(form.F, aes(x = SL_log, y = HW_log)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
F.reg.form.HW.plot
## `geom_smooth()` using formula = 'y ~ x'

F.reg.form.SnL <- lm(form.F$SnL_log ~ form.F$SL_log)
F.sd.form.SnL <- rstandard(F.reg.form.SnL)
F.reg.form.SnL.plot <- ggplot(form.F, aes(x = SL_log, y = SnL_log)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
F.reg.form.SnL.plot
## `geom_smooth()` using formula = 'y ~ x'

F.reg.form.OL <- lm(form.F$OL_log ~ form.F$SL_log)
F.sd.form.OL <- rstandard(F.reg.form.OL)
F.reg.form.OL.plot <- ggplot(form.F, aes(x = SL_log, y = OL_log)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
F.reg.form.OL.plot
## `geom_smooth()` using formula = 'y ~ x'

F.reg.form.FLA <- lm(form.F$FLA ~ form.F$SL)
F.sd.form.FLA <- rstandard(F.reg.form.FLA)
F.reg.form.FLA.plot <- ggplot(form.F, aes(x = SL, y = FLA)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
F.reg.form.FLA.plot
## `geom_smooth()` using formula = 'y ~ x'

##### MEX #####
F.reg.mex.D <- lm(mex.F$D ~ mex.F$SL)
F.sd.mex.D <- rstandard(F.reg.mex.D)
F.reg.mex.D.plot <- ggplot(mex.F, aes(x = SL, y = D)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
F.reg.mex.D.plot
## `geom_smooth()` using formula = 'y ~ x'

F.reg.mex.P1 <- lm(mex.F$P1 ~ mex.F$SL)
F.sd.mex.P1 <- rstandard(F.reg.mex.P1)
F.reg.mex.P1.plot <- ggplot(mex.F, aes(x = SL, y = P1)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
F.reg.mex.P1.plot
## `geom_smooth()` using formula = 'y ~ x'

F.reg.mex.P2.L <- lm(mex.F$P2.L ~ mex.F$SL)
F.sd.mex.P2.L <- rstandard(F.reg.mex.P2.L)
F.reg.mex.P2.L.plot <- ggplot(mex.F, aes(x = SL, y = P2.L)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
F.reg.mex.P2.L.plot
## `geom_smooth()` using formula = 'y ~ x'

F.reg.mex.P2.R <- lm(mex.F$P2.R ~ mex.F$SL)
F.sd.mex.P2.R <- rstandard(F.reg.mex.P2.R)
F.reg.mex.P2.R.plot <- ggplot(mex.F, aes(x = SL, y = P2.R)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
F.reg.mex.P2.R.plot
## `geom_smooth()` using formula = 'y ~ x'

F.reg.mex.A <- lm(mex.F$A ~ mex.F$SL)
F.sd.mex.A <- rstandard(F.reg.mex.A)
F.reg.mex.A.plot <- ggplot(mex.F, aes(x = SL, y = A)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
F.reg.mex.A.plot
## `geom_smooth()` using formula = 'y ~ x'

F.reg.mex.P1.R <- lm(mex.F$P1.R ~ mex.F$SL)
F.sd.mex.P1.R <- rstandard(F.reg.mex.P1.R)
F.reg.mex.P1.R.plot <- ggplot(mex.F, aes(x = SL, y = P1.R)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
F.reg.mex.P1.R.plot
## `geom_smooth()` using formula = 'y ~ x'

F.reg.mex.LLSC <- lm(mex.F$LLSC ~ mex.F$SL)
F.sd.mex.LLSC <- rstandard(F.reg.mex.LLSC)
F.reg.mex.LLSC.plot <- ggplot(mex.F, aes(x = SL, y = LLSC)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
F.reg.mex.LLSC.plot
## `geom_smooth()` using formula = 'y ~ x'

F.reg.mex.SALL <- lm(mex.F$SALL ~ mex.F$SL)
F.sd.mex.SALL <- rstandard(F.reg.mex.SALL)
F.reg.mex.SALL.plot <- ggplot(mex.F, aes(x = SL, y = SALL)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
F.reg.mex.SALL.plot
## `geom_smooth()` using formula = 'y ~ x'

F.reg.mex.SBLL <- lm(mex.F$SBLL ~ mex.F$SL)
F.sd.mex.SBLL <- rstandard(F.reg.mex.SBLL)
F.reg.mex.SBLL.plot <- ggplot(mex.F, aes(x = SL, y = SBLL)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
F.reg.mex.SBLL.plot
## `geom_smooth()` using formula = 'y ~ x'

F.reg.mex.SBDF <- lm(mex.F$SBDF ~ mex.F$SL)
F.sd.mex.SBDF <- rstandard(F.reg.mex.SBDF)
F.reg.mex.SBDF.plot <- ggplot(mex.F, aes(x = SL, y = SBDF)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
F.reg.mex.SBDF.plot
## `geom_smooth()` using formula = 'y ~ x'

F.reg.mex.BD <- lm(mex.F$BD_log ~ mex.F$SL_log)
F.sd.mex.BD <- rstandard(F.reg.mex.BD)
F.reg.mex.BD.plot <- ggplot(mex.F, aes(x = SL_log, y = BD_log)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
F.reg.mex.BD.plot
## `geom_smooth()` using formula = 'y ~ x'

F.reg.mex.CPD <- lm(mex.F$CPD_log ~ mex.F$SL_log)
F.sd.mex.CPD <- rstandard(F.reg.mex.CPD)
F.reg.mex.CPD.plot <- ggplot(mex.F, aes(x = SL_log, y = CPD_log)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
F.reg.mex.CPD.plot
## `geom_smooth()` using formula = 'y ~ x'

F.reg.mex.CPL <- lm(mex.F$CPL_log ~ mex.F$SL_log)
F.sd.mex.CPL <- rstandard(F.reg.mex.CPL)
F.reg.mex.CPL.plot <- ggplot(mex.F, aes(x = SL_log, y = CPL_log)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
F.reg.mex.CPL.plot
## `geom_smooth()` using formula = 'y ~ x'

F.reg.mex.PreDL <- lm(mex.F$PreDL_log ~ mex.F$SL_log)
F.sd.mex.PreDL <- rstandard(F.reg.mex.PreDL)
F.reg.mex.PreDL.plot <- ggplot(mex.F, aes(x = SL_log, y = PreDL_log)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
F.reg.mex.PreDL.plot
## `geom_smooth()` using formula = 'y ~ x'

F.reg.mex.DbL <- lm(mex.F$DbL_log ~ mex.F$SL_log)
F.sd.mex.DbL <- rstandard(F.reg.mex.DbL)
F.reg.mex.DbL.plot <- ggplot(mex.F, aes(x = SL_log, y = DbL_log)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
F.reg.mex.DbL.plot
## `geom_smooth()` using formula = 'y ~ x'

F.reg.mex.HL <- lm(mex.F$HL_log ~ mex.F$SL_log)
F.sd.mex.HL <- rstandard(F.reg.mex.HL)
F.reg.mex.HL.plot <- ggplot(mex.F, aes(x = SL_log, y = HL_log)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
F.reg.mex.HL.plot
## `geom_smooth()` using formula = 'y ~ x'

F.reg.mex.HD <- lm(mex.F$HD_log ~ mex.F$SL_log)
F.sd.mex.HD <- rstandard(F.reg.mex.HD)
F.reg.mex.HD.plot <- ggplot(mex.F, aes(x = SL_log, y = HD_log)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
F.reg.mex.HD.plot
## `geom_smooth()` using formula = 'y ~ x'

F.reg.mex.HW <- lm(mex.F$HW_log ~ mex.F$SL_log)
F.sd.mex.HW <- rstandard(F.reg.mex.HW)
F.reg.mex.HW.plot <- ggplot(mex.F, aes(x = SL_log, y = HW_log)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
F.reg.mex.HW.plot
## `geom_smooth()` using formula = 'y ~ x'

F.reg.mex.SnL <- lm(mex.F$SnL_log ~ mex.F$SL_log)
F.sd.mex.SnL <- rstandard(F.reg.mex.SnL)
F.reg.mex.SnL.plot <- ggplot(mex.F, aes(x = SL_log, y = SnL_log)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
F.reg.mex.SnL.plot
## `geom_smooth()` using formula = 'y ~ x'

F.reg.mex.OL <- lm(mex.F$OL_log ~ mex.F$SL_log)
F.sd.mex.OL <- rstandard(F.reg.mex.OL)
F.reg.mex.OL.plot <- ggplot(mex.F, aes(x = SL_log, y = OL_log)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
F.reg.mex.OL.plot
## `geom_smooth()` using formula = 'y ~ x'

F.reg.mex.FLA <- lm(mex.F$FLA ~ mex.F$SL)
F.sd.mex.FLA <- rstandard(F.reg.mex.FLA)
F.reg.mex.FLA.plot <- ggplot(mex.F, aes(x = SL, y = FLA)) + 
  geom_point() +
  stat_smooth(method = "lm", col = "red") +
  stat_cor(label.y = 10)
F.reg.mex.FLA.plot
## `geom_smooth()` using formula = 'y ~ x'